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Preface 


This  thesis  topic  was  proposed  by  Air  Force  Wright 
Aeronautical  Laboratories,  WPAFB,  Ohio.  It  is  an  analytic 
study  of  the  effects  of  surface  roughness  on  compressible 
turbulent  boundary  layer.  It  interested  me  because  of  the 
mysterious  nature  of  the  turbulence,  and  I  decided  that  I 
ought  to  undertake  this  investigation  and  attempt  to  learn 
as  much  as  I  could  about  this  subject.  The  topic  is  also 
of  current  interest  to  USAF. 

I  take  this  opportunity  to  thank  "Almighty  Allah" 
who  gave  me  the  strength  to  accomplish  this  project.  I 
wish  to  express  my  sincere  gratitude  io  my  thesis  advisor, 
Dr.  Urmilla  Ghia,  whose  understanding  oi  the  subject  was 
very  impressive.  Her  continued  guidance,  contagious 
enthusiasm  and  above  all  unlimited  patience  throughout 
all  phases  of  this  study  were  invaluable.  I  also  wish 
to  thank  Dr.  William  Hankey,  Jr. ,  for  his  sponsorship  and 
continued  assistance-.  I  could  not  give  this  paper  a 
proper  preface  without  acknowledging  Capt  James  K.  Hodge 
for  his  constant  interest.  His  support,  be  that  in  the 
form  of  enlightening  criticism  or  compliments,  often  came 
as  insight  during  times  of  frustrations.  I  will  always 
remember  with  fondness  his  friendly  and  calm  nature,  and 
his  devotion  to  duty.  I  also  appreciate  the  diligence  of 
Mrs.  Cindy  Boone  who  had  immense  patience  to  type  this 
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manuscript  flooded  with  mathematical  equations  and  greek 

symbols . 

Last  but  not  the  least,  I  express  individual  thanks 
to  my  sweet  wife,  for  her  loving  support,  assistance  and 
understanding  with  which  she  took  the  late  hours  and  the 
weekends  not  only  during  the  thesis,  but  during  my  whole 
stay  at  AFIT.  I  also  thank  my  children,  Hannah,  Ahad 
and  Fahad  for  their  patience  with  !'Baba"  as  he  reached 
for  his  goal  of  achieving  higher  learning. 
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Abstract 


This  study  followed  the  work  of  Dr.  Anthony  Fiore,  of 

Air  Force  Wright  Aeronautical  Laboratories,  Wright-Patterson 

Air  Force  Base,  Ohio.  Dr.  Fiore  had  carried  out  an  expert¬ 
ly  fr'viPV, -  rovfl 

\  raentalr study  ,i5T  the  effect  of  surface  roughness  on  the 

1  Q ? 

turbulent  boundary  layery^ (£dx-&rrf&  Fortran  code,  ITRACT,  ^ 
written  primarily  by  Dr.  Shang,  that  solved  for  the  charac¬ 
teristics  of  a  laminar,  transitional  and  turbulent  boundary 
layer  on  smooth  surfaces.  ^Fke-.puiqoo-se  -  of1'  the  present 
study  4as  -to  "investigate5’the  influence  of  surface  roughness 
on  a  compressible  turbulent  boundary  layer  and  then 
extendijthe  usefulness  of  hhe  existing  computer  code^  ITRACT,  ^ 
by  including  in  it  the  optional  capability  of  rough-surface 
boundary-layer  calculations.  ^ 

To  achieve  this  objective,  the" ^surface  roughness  was 
represented  by  distributed  sources  and  sinks  in  the  appro¬ 
priate  governing  equations.  The  most  important  term  is  a 
sink  term  in  the  mean  momentum  equation,  representing 

9 

the  form  drag  due  to  the  roughness  element.  T-he  governing 

' 

boundary-layer  equations  for  continuity,  momentum,  and 


energy  were  derived  in  a  form  to  account  for  t-hre  blockage 
effect^  jdue  to  thoroughness  elements.  The  modified  govern- 
ing  equations  were  then  transformed  using  transformation 

r 0,1  s II'  5  , 

Probs tein-Elliott  and  Levy-Lees^  The  resulting  equa¬ 
tions,  with  appropriate  boundary  conditions,  were  solved 

] 

by  finite-difference  techniques  to  determine  the  non-  -> 


dimensional  velocity  components  and  temperature  at  a 

ivV'’ 

finite  number  of  nodes  in  the  boundary-layer f field . of  the 

■  A\ 

flow. 

To  establish  the  authenticity  of  the  original  code, 
ITRACT,  some  smooth-surface  results  were  first  computed 
and  compared  with  the  experimental  data  of  Dr.  Fiore  and 
Dr.  Cole  for  turbulent  flows  over  smooth  surfaces.  .  After 
the  accuracy  of  the  code  was  established  for  smooth-surfa 
calculations,  the  code  was  modified  in  order  to  render  it 
capable  of  predicting  the  influence  of  surface  roughness 
on  compressible  turbulent  flow.  The  modified  code  was 
then  used  to  obtain  results  for  rough-surface  boundary 
layers  and  the  computed  results  were  compared  with  the 
experimental  data  of  Dr.  Fiore  for  the  case  of  supersonic 
flow  over  a  rough  flat  plate.  The  agreement  between  the 
computed  and  the  measured  velocity  profiles  was  quite 
satisfactory.  The  corresponding  temperature  profiles 
agreed  well  everywhere,  except  very  near  the  wall;  a 
possible  reason  for  this  discrepancy  is  offered  in  this 
study.  Unlike  previous  studies  of  rough-surface  boundary 
layers,  the  present  study  makes  no  modification  to  the 
turbulence  model  employed. 
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AN  ANALYTIC  STUDY  OF  THE  EFFECT  OF  SURFACE  ROUGHNESS 
ON  A  COMPRESSIBLE  TURBULENT  BOUNDARY  LAYER 

I .  Introduction 

Preliminaries  and  Problem  Analysis 

Surface  roughness  can  play  an  important  role  in 
turbulent  boundary-layer  skin  friction  and  heat  transfer 
for  many  high-speed  flight  applications.  An  understand- 
inf  of  the  roughness  effect  is  essential  for  accurate  de¬ 
sign  prediction  in  a  wide  variety  of  applications,  includ¬ 
ing  ships,  aircraft,  compressor  blades,  turbine  blades, 
missiles  and  re-entry  vehicles.  For  example,  the  NASA  Space 
Shuttle  Program  studied  roughness  as  it  augments  heating. 

At  low  flight  altitude,  the  thickness  of  boundary  layer  on 
the  blunted  nose  region  of  ypersonic  re-entry  vehicle  can 
easily  be  less  than  the  inherent  surface  roughness  of  prac¬ 
tical  heat  shield  materials,  and  roughness  dominates  the 
heat  transfer  characteristics.  Also  recent  experiments 
(Ref.  1)  have  shown  that  surface  roughness  alone  can  sig¬ 
nificantly  influence  the  control  effectiveness  of  maneu¬ 
vering  vehicle.  Numerical  analysis  and  computer  programs 
have  been  developed  over  the  last  twenty  years  to  solve 
turbulent  boundary  layers  over  smooth  surfaces.  A  defi¬ 
ciency,  however,  exists  with  these  programs  as  regards 
inclusion  of  roughness  effects.  Data  has  been  accumulated 
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over  the  last  few  years  in  this  area  but  a  need  for  a  com¬ 
puter  code  to  accurately  match  this  data  still  exists. 

The  purpose  of  this  study  was  to  investigate  the  influence 
of  surface  roughness  on  a  compressible  turbulent  boundary 
layer  and  then  to  modify  an  existing  computer  code  for 
turbulent  boundary  layer  over  smooth  surface  to  account 
for  roughness  effects. 

Literature  Survey 

Most  available  models  for  analyzing  the  influence  of 
surface  roughness  on  boundary-layer  behavior  are  essentially 
extension  of  Nikuradse's  study  (Ref.  2)  of  pipes  roughened 
with  sand  and  application  of  Nikuradse's  results  to  flat 
plates  by  Prandtl  and  Schlichting  (Ref.  3).  Several  cor¬ 
relations  have  been  proposed  to  relate  real  surface- 
roughness  heights,  spacing  and  geometries  to  an  equivalent 
sand-grain  roughness  height  so  that  Nikuradse's  data  can 
be  used.  Examples  of  such  correlations  can  be  found  in 
White  and  Grabow  (Ref.  4 s 1 53-164).  Dvorak  (Ref.  5:1752- 
1759)  used  integral  methods  in  which  the  skin-friction 
coefficient  was  specified  as  a  function  of  boundary-layer 
thickness  and  roughness  height.  Using  this  specification, 
the  moment  equations  were  solved  for  the  momentum  and  dis¬ 
placement  thickness.  Chen  (Ref.  6:623-629)  extended  this 
approach  to  predict  heat  transfer,  by  using  a  Stanton  num¬ 
ber  correlation  derived  from  the  subsonic  data  by  Owean 
and  Thomson  (Ref.  7:321-334).  In  this  approach,  the  stag- 
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nation  enthalpy  profile  was  assumed  to  have  the  same  shape 
as  the  velocity  profiles.  A  similar  model  has  recently 
been  developed  for  re-entry  vehicles  by  Dahm  et  al.  (Ref.  8). 
Here  again  a  momentum  integral  approach  is  used,  with  the 
skin-friction  and  heat-transfer  coefficients  based  on  cor¬ 
relations  of  the  low-speed  data  of  Healzer  et  al.  (Ref.  9) 
and  flat  plate  measurements  of  Reda  (Ref.  10)  at  a  macb 
number  of  2.9.  The  roughness  augmentation  of  heat  transfer 
was  found  to  be  about  60  percent  of  the  skin-friction  aug¬ 
mentation.  More  recently,  effects  of  surface  roughness 
have i been  evaluated  by  differential  methods.  Cebeci  and 
Chang  (Ref  11:730-735)  numerically  solved  the  incompies- 
sible  boundary  layer  equations  employing  an  algebraic 
eddy  viscosity  formulation  modified  for  surface  rough¬ 
ness.  The  modification  was  based  on  Rotta's  (Ref.  12:1- 
219)  model,  which  displaces  the  normal  coordinate  of  the 
rough-wall  velocity  profile.  An  expression  for  this  dis¬ 
placement  and  the  resulting  mixing  length  is  given  by 
Cebeci  and  Chang  (Ref.  11)  as  a  function  of  an  equivalent 
sand-grain  roughness  height.  Emphasizing  compressible 
flows  for  a  variety  of  edge  and  wall  conditions,  Hodge 
and  Adams  (Ref.  13)  numerically  solved  the  flow  equations 
together  with  the  equation  of  kinetic  energy  of  turbulence. 
Roughness  effects  were  accounted  for  by  inclusion  of  a 
form  drag  term  in  the  momentum  equation  and  by  modification, 
based  on  the  results  of  Healzer,  et  al.,  of  several  of  the 
nine  empirical  constants  in  the  turbulence  model.  A 
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somewhat  more  involved  approach  was  taken  by  Saffman  and 
Wilcox  (Ref.  142541-546)  utilizing  a  two-equation  turbulence 
model.  However,  the  effect  of  roughness  was  treated  rather 
empirically  by  making  the  boundary  conditions  for  the  pseudo- 
vorticity  at  the  wall  a  function  of  the  roughness  height. 

This  dependence  was  derived  so  as  to  fit  the  observed  vari¬ 
ation  of  the  "law  of  the  wall"  velocity  deficit  with  rough¬ 
ness.  Some  encouraging  profiles  were  also  computed  for  the 
mean  and  fluctuating  velocities.  However,  heat  transfer 
was  again  determined  by  invoking  a  Reynolds  analogy  with 
the  skin  friction.  Finson  and  Clark  (Ref.  15:3-6)  pre¬ 
sented  a  technique  which  accounted  fo^  the  surface  rough¬ 
ness  by  calculating  the  form  drag  contribution  by  individ¬ 
ual  elements.  Following  the  approach  of  Finson  and  Clark, 
Christoph  ana  Fletcher  (Ref.  16:509-510)  used  a  two-layer 
algebraic  mixing  length  mod<=!  that  explicitly  accounts  for 
mass  addition  and  surface  roughness,  in  addition  to  the 
modification  of  the  boundary-layer  equations  as  suggested 
by  Finson  and  Clark. 


Scope  of  Present  Stud1 


The  Flight  Dynamics  Laboratory  possessed  a  digital 
computer  code  called  ITRACT,  which  computed  the  character¬ 
istics  of  laminar  and  turbulent  boundary  layers  for  either 
planar  or  axisymmetric  flow  over  smooth  surfaces.  The 
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purpose  of  the  present  study  was  to  modify  this  code  for 
inclusion  of  the  effects  of  surface  roughness  on  compres¬ 
sible  flow.  Historically,  roughness  effects  have  been 
modelled  by  a  law-of-the-wall  velocity  profile  expression 
in  terms  of  an  equivalent  sand-grain  roughness  height. 
Physically,  the  equivalent  sand-grain  roughness  concept 
is  not  very  satisfying  since  an  equivalent  sand-grain 
roughness  height  must  be  contrived  for  real  roughness 
heights,  spacing  and  geometries.  A  physically  more 
meaningful  method  is  that  employed  'y  Finson  an i  Clark 
and  followed  by  Christoph  and  Pletcher.  In  the  present 
study,  the  same  approach  is  followed  but  without  invoking 
any  modificatjon  of  the  turbulence  model.  Roughness  is 
represented  by  distributed  sources  and  sinks  in  one 
appropriate  governing  equations.  The  most,  important 
term  is  a  sink  term  in  the  mean  momentum  equaoirr.  rep¬ 
resenting  the  form  drag  on  the  roughness  e1  sm.er.o  s ,  The 
governing  boundary-layer  equations  are  cast  in  a  form  to 
account  for  the  blockage  effect  due  to  the  roughness 
elements.  Accordingly,  the  fluxes  along  the  streamwise 
direction  are  multiplied  by  (1  -  D(y)/£)  where  I(y)  is 
the  element  diameter  at  height  y  and  Z  the  average  centre 
to  centre  spacing  of  the  element  (Fig  1).  Fluxes  in 
a  direction  normal  to  the  streamwise  direction  are  multi- 
plied  by  (l  -  ttD  /4&  ).  Here,  the  shape  of  the  elements 
has  been  restricted  to  circular  cross-sections  cnly,  but 
the  modified  code  would  eventually  have  the  provision  to 


predict  the  roughness  effects  due  to  rectangular  elements 
as  well.  The  computed  results  obtained  by  employing  a 
Reynolds  stress  turbulence  model  in  combination  with  a 
drag  description  for  the  effect  of  the  roughness  elements 
on  the  flow  is  compared  against  relevant  data  to  establish 
the  validity  and  accuracy  of  the  theory  and  to  offer  ex¬ 
planations  for  the  observed  trends. 

The  model  used  in  this  study  is  aimed  entirely  at 
distributed  roughness,  i.e.,  three-dimensional  roughness, 
appropriate  to  the  vast  majority  of  practical  applications. 
Two-dimensional  roughness  such  as  machined  grooves  normal 
to  the  flow  direction  are  not  considered  here.  The  two 
types  of  roughness  may  yield  qualitatively  similar  trends 
in  terms  of  roughness  height,  spacing,  etc.,  but  a  sub¬ 
stantial  difference  in  the  nature  of  the  flow  may  be 
expected.  This  model  makes  the  basic  assumption  that  the 
forces  on  roughness  elements  can  be  viewed  as  form  drag. 
This  implicitly  requires  that  the  flow  approaching  an 
individual  element  be  attached,  whereas  ivith  2-D  rough¬ 
ness  elements,  cavity  flow  is  likely  to  prevail  immediately 
downstream  of  the  roughness  elements. 

The  next  step  in  the  study  was  to  learn  as  much  about 
the  computer  code,  ITRACT,  as  possible.  This  step  included 
a  study  of  the  equations  of  motion,  continuity,  and  energy, 
together  with  the  perfect  gas  law  and  Sutherland's  viscos¬ 
ity  law  needed  for  a  boundary-layer  calculation.  It  also 
included  the  comparison  of  the  computer  code  results  for 
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turbulent  flow  over  smooth  surfaces  to  experimental  data 
to  establish  the  authenticity  of  the  code. 

After  the  accuracy  of  tae  present  code  was  established, 
the  next  step  was  to  modify  the  coda  in  order  to  render  it 
capable  of  predicting  the  effect  of  surface  roughness  on 
compressible  turbulent  flow.  The  modified  code  was  verified 
by  comparing  the  results  predicted  by  modified  computer  code 
with  the  experimental  data  available  at  the  Flight  Dynamics 
Laboratory  for  the  case  of  supersonic  flow  over  rough  flat 
plate.  This  comparison  will  also  establish  the  extent  to 
which  the  modification  o'f  the  turbulence  model  is,  or  is 
not,  needed  for  accurate  prediction  of  rough-surface 
boundary  layer  flows. 

The  major  accomplishment  of  this  study  is  the  exten¬ 
sion  of  the  usefulness  of  the  existing  computer  code, 

ITRACT,  by  including  in  it  the  optional  capability  of 
rough-surface  boundary-layer  calculations. 


II .  Analysis  of  Problem 
Governing  Equations 

This  section  presents  the  governing  equations  for  the 
compressible  turbulent  boundary  layer  together  with  the 
required  boundary  conditions.  In  their  final  form,  govern¬ 
ing  equations  include  the  effect  of  surface  roughness. 

The  eddy  viscosity  and  eddy  conductivity  models  used  to 
represent  the  apparent  turbulent  shear  and  heat  flux  terms 
appearing  in  the  mean-flow  boundary-layer  equations  are 
discussed  in  the  latter  part  of  this  section. 

Coordinate  System 

The  coordinate  system  employed  is  shown  in  Figure  2. 

The  boundary-layer  coordinate  system  is  denoted  by  x  and  y 
which  are  tangent  and  normal  to  the  surface,  respectively. 
The  origin  of  the  boundary-layer  coordinate  system  x,  y, 
as  well  as  that  of  the  body  coordinate  system  z,  r,  for 
axisymmetric  configurations  is  located  at  the  stagnation 
point  for  blunt  bodies,  and  at  the  leading  edge  for  sharp- 
tipped  bodies.  The  velocity  components  u  and  v  are  oriented 
along  the  x  and  y  directions,  respectively .  Transverse 
curvature  terms  are  retained  because  of  their  importance 
in  the  development  of  boundary-layer  flow  over  slender 
bodies  of  revolution  where  the  boundary-layer  thickness 
may  become  of  the  order  of  the  body  radius  r  .  The  angle 
<t>  is  the  angle  between  z  axis  and  local  tangent  evaluated 
at  (x, 0 ) . 
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Differential  Equations 

The  flow  of  a  compressible,  viscous,  heat  conducting 
fluid  is  mathematically  described  by  the  continuity  equa¬ 
tion,  the  Navier-Stokes  equations  and  the  energy  equation, 
together  with  an  equation  of  state,  a  heat  conductivity  law 
and  a  viscosity  law.  For  flows  at  large  Reynolds  number, 
Prandtl  (Ref.  17)  has  shown  that  the  Navier-Stokes  equations 
and  the  energy  equation  can  be  simplified  to  a  form  now 
recognized  as  the  compressible  boundary-layer  equations^ 
These  equations  may  be  written  as  follows: 

Continuity : 

9  (r^pu)  +  (rJ'pv)  =  0  (1 

9x  9y 

Momentum : 


[’ll  3u  + 

v  3u”| 

=  -£e  + 

1  3  I"  r'1 

fv  lu\l 

[  9x 

ayj 

dx 

r.i  sy[ 

Energy: 

fu  9  (CpT)  +  v  9  (Cpt)  1 
p  l  3x  3y  ,| 


u  do 
dx 


+  1_  9_  |"rJ  IU  9_  ( CpT )  ]  +  p/9u\ 
rj  9y  Cp  9y  J  \9y  ) 


Osborne  Reynolds,  who  was  first  to  observe  and  study 
the  phenomenon  of  transition  from  laminar  to  turbulent  flow, 
assumed  that  the  instantaneous  fluid  velocity  satisfied  the 


Navier-Stokes  equations  and  that  the  instantaneous  velocity 
(u^)  could  be  considered  to  consist  of  a  mean  (time  averaged) 
component  u  and  a  fluctuating  component  u',  i.e., 

ui(xi,t)  =  C )  +  u1"(xi,t)  where  i  =  1,2,3  (4) 

In  order  to  obtain  the  conservation  equations,  the  instan¬ 
taneous  quantities  in  the  equations  (1)  to  (3)  were  re¬ 
placed  by  their  mean  and  their  fluctuating  quantities. 

By  taking  the  time  average  of  the  various  terms  appearing 
in  these  equations  and  making  the  boundary-layer  assump¬ 
tions,  the  following  mean  continuity,  mean  momentum  and 
mean  energy  equations  were  derived  (Ref.  18,  145,  216). 

9  (r^pu)  +  3__  f  rJ* p  (  v  +  p  ’v  '  0=0  ( 5 ) 

3x  3y  L  p  'J 


Momentum : 


These  equations  are  identical  to  those  for  laminar 
flows,  with  the  exception  of  the  correlations  of  the  turbu¬ 
lent  fluctuating  quantities  which  represent  the  apparent 
mass,  shear,  and  heat-flux  terms  caused  by  turbulence. 

These  fluctuating  quantities  were  incorporated  throu  gb 
mathematical  modelling.  The  apparent  mass  flux  term  p 'v  * , 
the  apparent  shear  stress  term  pu >v J  and  the  apparent  heat 
flux  term  Cppv '  are  represented  by  a  new  velocity  compo¬ 
nent  v,  an  eddy  viscosity  e,  and  an  eddy  conductivity  Km. 
respectively. 

These  terms  were  defined  by  the  following  relationships: 


The  turbulent  Prandtl  number  is  expressed,  in  a  manner 
analogous  to  the  laminar  Prandtl  number  expressed  in  terms 
of  viscosity  and  eddy  conductivity  as: 


(11  ) 


To  this  set  of  equations  the  following  perfect  gas  law 


and  the  viscosity  relation  of  Sutherland  were  also  added. 


r.’V  ■■-_  ^  «r^  ■r*.  a  *r;  «r\  *-.  tj  *"i  WT  V 


.s  v,"  '.^v^  ;^» v* '."» 


.»*;*•'* 


Perfect-gas  law: 


P  = 


Cp 


pT 


Viscosity  Law: 


(air  only) 


(12) 


(13) 


iV) 

s*> 


where  p  denotes  the  viscosity  at  the  reference  T  and  ’S' 

@ 

is  a  constant.^  This  relation  is  approximated  in  theoretical 
calculation  by  the  simpler  power-law: 

M_  =  /  T\W  0 . 5<w<1 . 0  (1 

ue  \  e  / 


It  has  been  found  that  Sutherland's  formula  cnn  be  approxi¬ 
mated  at  high  temperature  by  adopting  values  of  w  between 
0.5  and  0.75,  whereas  at  lower  temperature  the  value  of  io=1.0 
appears  to  be  adequate. 

Rough-Wall  Boundary  Layer  Model 

The  basic  model  for  the  rough-wall  boundary-layer  is 
the  same  as  described  thus  far  for  smooth  wall  boundary- 
layer  flow  and  free-shear  flows.  The  roughness  model  makes 
the  basic  assumption  that  the  additional  forces  due  to  the 
roughness  elements  can  be  viewed  as  form  drag.  This  im¬ 
plicitly  requires  that  the  flow  approaching  an  individual 
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element  be  attached.  This  model  is  more  appropriate  for 
distributed  roughness,  i.e.,  3-D  roughness  elements,  since 
cavity  flow  is  likely  to  prevail  with  2-D  roughness. 

The  rough  surface  is  idealized  as  being  made  up  of 
identical  elements.  The  bottom  of  the  elements,  or  the 
underlying  smooth  wall,  is  at  y=0  (Fig  1).  The  element 
height  is  K,  and  l  is  the  average  element  spacing  measured 
from  centre  to  centre  of  two  adjacent  elements.  The  total 
number  of  elements  per  unit  area  is  given  byl/il  .  The 
analysis  presented  here  is  for  a  case  of  roughness  elements 
with  circular  cross-section  at  all  heights,  with  D(y)  de¬ 
noting  the  diameter  of  the  element  at  height  y  for  0<y<K, 
but  any  general  shape  may  be  specified.  Then,  viewing  flow 
around  the  element  at  height  y  sc  two  dimensional,  the  form 

drag  between  y  -  6y  and  y  +  5y_  is: 

2  2 

1  pu2  CD  D(y)6y 

2 

where  is  the  form  drag  coefficient.  To  relate  this  to 
drag  per  unit  volume,  it  is  noted  that  there  are  Z~  ele¬ 
ments  per  unit  area,  so  that  the  appropriate  differential 
2 

volume  is  l  6y  and,  therefore,  the  sink  term  for  mean 
momentum  is: 

Ru  =  -1  pu2  CD  M  ( 

jr 

The  drag  coefficient  could  be  specified  to  be  equal  to  unity 
(Cp  =  1),  appropriate  to  infinite  circular  cylinders  (two- 
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dimensional  elements)  at  local  Reynolds  number  above  the 
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Stokes-flow  regime.  However,  lower  values  such  as  =  0.6 
are  more  appropriate  for  finite  elements  (3-D  elements) 


such  as  cones,  hemisphere,  etc. 

In  addition,  there  should  be  source  terms  in  the  equa¬ 
tion  governing  turbulent  kinetic  energy  and  its  dissipation 
rate  in  order  to  describe  the  tendency  of  roughness  to  in¬ 
crease  the  velocity  fluctuation.  Specification  of  these 
terms  is  more  speculative;  their  contributions  are  generally 
smaller  than  the  natural  production  terms  and,  hence,  less 


critical  to  the  model.  These  terms  are  not  very  important 
compared  to  the  indirect  effect  of  roughness  to  increase 
the  turbulent  energy  by  increasing  the  mean  shear.  Except 
in  the  Stokes  flow  regime,  heat  transfer  to  an  element  should 
be  small.  Therefore,  the  only  roughness  term  appearing  in 
the  thermal  equations  is  a  source  term  in  the  mean  static 
enthalpy  equation.  This  term  is  constructed  such  that,  in 
combination  with  the  sink  term  (Equation  15)  for  drag,  the 
total  enthalpy  is  not  altered.  Accordingly,  the  mean  static- 
enthalpy  differential  equation  must  include  a  source  term  R^, 
defined  as: 

Rh  =  pu3  CD  V(y)/i2  ( 

The  detailed  derivation  of  this  term  is  given  in  Appendix  A. 

If  no  further  modification  is  made  in  the  governing 
equations,  it  is  implied  that  the  roughness  elements  are 
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assumed  to  occupy  no  space  (Ref.  25:2).  This  assumption 
becomes  progressively  worse  as  the  roughness  density  in¬ 
creases,  therefore,  the  model  has  been  extended  to  account 
for  the  blockage  effects  due  to  the  roughness  elements. 
Accordingly,  the  boundary-layer  equations  are  derived  in  a 
form  to  account  for  the  blockage  effect.  This  is  done  in 
the  following  manner: 

(a)  At  a  given  height  y,.the  fraction  of  flow  area  in 
the  x-direction,  that  is  open  to  the  flow,  is  { 1 -D(y)/£} , 
hence,  fluxes  in  the  streamwise  direction,  represented  via 
the  convective  operator  pu  3/3x,  are  multiplied  by  this 
factor. 

(b)  Fluxes  across  a  surface  area  whose  normal  is  in 
the  y-direction,  should  be  multiplied  by  { 1  -  D2(y)/4 12}. 
However,  the  roughness  terms  discussed  above  are  already 
based  on  the  total  volume,  rather  than  the  available  flow 
volume,  and  need  no  such  factor. 

With  the  incorporation  of  these  modifications,  the 
conservation  equations  (5),  (6),  and  (7)  are  recast  as 
follows : 

Continuity : 

f(y)  JL_  (r^'pu)  +  d_  (r^pv)  =  0 
3x  9y 


Momentum : 


viv  s,' v  v.v%n  .  v  V.  V  *.  •/  VX  •.  \  v* 


w  7?  ’  t-.  w-"  i  t"  irr  t*"'  ■ 


.  r-. 
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Energy: 


+  1  1 
r5  BlyT 

3_  jB(y)rJ  Ju 

3u  -  pu'vH  I 

3y  JJ 

-1  pu2  CD  DW 

0  O 

(18) 

2  B  (y )  S,2 

3 ( CpT )  +  v 

3  ( CpT )  =  f  ( ./ )  u 

d£ 

3x 

3y 

dx 

© 


BT?T: 


3 

B(y)rj 

[KJl  3(CpT)i| 

+  U  f  3u  T 

3y 

- 

LCp  ay  JJ 

L  3y  J 

3_ 

B  ( y )  r  ^ 

(-Cppv’T')l 

-  pu'v'  f3u-| 

9y 

J 

LsyJ 

+  1_  pu3  CD  D(y) 
2  B ( y ) £2 


(19) 


In  equations  (17-19): 


B(y)  =  {1  -  7rD2(y)/4il2 


(20) 


f(y)  =  (1  -  D(y)/Jl)/B(y) 


(21) 


>3 


Details  of  the  derivation  of  equations  (17)  through  (19) 
are  given  in  Appendix  B.  The  function  f(y)  contains  the 
main  effect  of  blockage.  If  a  stream  function  formulation 
were  incorporated,  f(y)  may  be  absorbed  in  the  definition  of 
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the  stream  function  that  is  introduced  to  eliminate  the 


normal  velocity. 


3fl  =  f(y)pu 
3y 


(22) 


M  =  -pv 
3x 


(23) 


However,  it  is  not  done  in  the  present  work.  It  should  be 
noted  that,  if  the  elements  are  packed  so  tightly  that  they 
are  touching  over  some  range  of  y,  then  3  =  l  and  f (y)  =  0 
over  that  range.  This  formulation  forces  the  velocity  to 
remain  zero  up  to  the  height  where  D<£  and  the  flow  is 
blocked.  In  such  cases,  Y  =  0  is  redefined  as  the  lowest 
point  where  the  flow  is  unblocked.  (Fig  1) 

A  major  advantage  of  this  model  is  that  solutions  are 
obtained  for  both  velocity  and  thermal  variables.  Heat 
transfer  is  obtained  directly  without  invoking  a  Reynolds 
analogy.  Finite-difference  solut-’  ns  are  obtained  using 
the  boundary  conditions  that,  (i)  Fluctuating  quantities 
are  zero  at  the  base  of  the  wall,  Y  =  0  and  (ii)  At  the 
outer  edge,  fluctuating  quantities  are  zero  and  the  flow 
variables  approach  the  f ree-streair.  values. 

Transformation  of  Boundary-Layer  Equations 

Equations  (17-19)#  which  are  expressed  in  the  surface- 
normal  coordinates  in  the  physical  plane,  require  starting 
profiles,  but  these  equations  are  singular  at  the  stagna- 


tion  point.  For  this  reason,  the  equations  are  transformed 
to  a  coordinate  system  that  removes  the  singularity  at  the 
stagnation  point,  stretches  the  coordinate  normal  to  the 
flow  direction,  thereby,  resulting  in  a  more  gradual  growtn 
of  the  boundary-layer  thickness  and  places  the  equations  in 
an  almost  two-dimensional  form  (Ref.  15).  A  combination  of 
the  Probstein  (Ref.  31 )  and  Levy-Less  was  used  in  this 
analytic  study.  The  transformed  coordinates  £  and  p  are 
defined  as : 

C(x)  =  0/xpe  ue  ue  r^dx  (24) 


and 


1 (x,y) 


u  r  2j 

e  o 
2C 


Q/ytJ’  dy 


(25) 


It  may  be  pointed  out  here  that  turbulent  boundary 
layers  are  characterized  by  two  length  scales,  namely,  the 
boundary-layer  thickness  and  the  wall-layer  thickness,  which 
are  quite  different  in  magnitude  and  vary  in  the  streamwiso 
direction  depending  upon  the  pressure  gradient,  wall  boundary 
conditions,  etc.,  thereby,  making  the  analysis  of  turbulent 
layers  more  complicated  than  the  laminar  boundary  layers 
where  generally  only  one  length  scale  is  present.  In  com¬ 
pressible  laminar  flow,  when  the  boundary-layer  equations 
are  expressed  in  terms  of  Levy-Lees  variables,  the  streamwise 
growth  of  the  boundary  layer  is  significantly  reduced  thereby 
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simplifying  the  numerical  solution  of  the  boundary  layer, 
riut  in  turbulent  flow  since  the  Levy-Lees  variables  do  not 
properly  capture  the  boundary  layer  thickness,  it  is  neces¬ 
sary  to  monitor  the  numerical  solution  and  add  computa¬ 
tional  points  in  the  outer  region  to  accommodate  the 
boundary-layer  growth.  To  overcome  this  problem,  it  would 
be  better  to  use  the  new  self-adaptive  coordinate  trans¬ 
formation  for  finite-difference  solution  of  turbulent 
boundary-layer  flows  presented  by  J.E.  Carter  et  al. 

tef.  32)  as  this  permits  a  uniform  mesh  to  be  used  in  the 
computational  coordinate  which  extends  across  the  layer. 
This  coordinate  transformation  uses  the  local  value  of  the 


skin  friction  to  scale  the  thickness  of  the  wall-layer 
region,  and  the  local  maximum  value  of  turbulent  viscosity 
tr  ale  the  boundary-layer  thickness. 

To  proceed  with  the  derivation  of  the  equations  used 
in  the  present  study,  the  dependent  variables  are  non- 
dimens.  inalized  according  to: 


F  =  u 


0  =  T_ 
T 


(26) 


Next,  the  relations  between  derivative  in  the  physical 
(x,y)  plane  and  in  the  transformed  plane  (£,n)  are  obtained 
as  follows : 


[fe] 


=  p  u  y  r 
K  e  e  e 


0  J  [  ft  ]  n  +  [fx]  [i]c 


(27) 


I 


and 


u  r 
e  o 


JtJ 


tt 


pA 


(28) 


In  Eqs  (27)  and  (28),  the  subscript  outside  the  sauare 
bracket  indicates  the  coordinate  that  is  held  constant  during 
the  indicated  differentiation  process. 

The  transformed  variable  for  the  normal  velocity  V  is 
obtained  as: 


V  = 


3L 


f ( y ) F  9a  +  roJtJpv 
9x  ZK 


p  u  u  r 
Ke  e  Me  o 


2J 


(29) 


The  detailed  derivation  of  equation  (29)  is  given  in 
Appendix  B. 

With  this,  the  final  working  form  of  the  governing 
system,  prior  to  linearization,  was  obtained  as  follows: 
(Ref.  Appendix  B) 

Continuity : 

Vn  t  f ( y ) F  +  2C  f(y)Fc  =  0  (30) 


Momen  ium : 


2?  f(y)FF?  +  0f(y)(F2-6)  +  VF^ 

=  _1  {B(y)t2jaF  }  -  $  (31) 

bTTT  n  n  ss 
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Energy : 


25  f (y)Fe.  +  ve; 


B(y)e 


l  9 
Ft  n 


] 


n 


+  t 


^JlaeFn^ 


$00aF' 

ss 


(32) 


Further,  the  additional  symbols  included  in  these  equations 
are  defined  as  follows: 
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Inner-Region  Model.1  The  eddy  viscosity  model  used  for 
the  inner  region  is  based  on  the  mixing-length  hypothesis 
of  Prandtl  (Ref.  21).  The  eddy  viscosity  for  this  region, 
referenced  to  the  molecular  viscosity,  may  be  expressed  as: 


(41  ) 


where  Z,  the  mixing  length,  may  be  written  as: 

I  =  K1  y 


(42) 


To  account  for  the  region  close  to  the  wall,  Van  Driest 
(Ref.  22)  suggested  a  modif ication  for  the  mixing  length  of 
Prandtl.  The  correct  form  for  the  mixing  length  in  the 
viscous  sub-layer  was  given  as: 

%  =  K.j y{  1  -  exp ( — y / A ) )  (43) 

where  the  exponential  term  is  due  to  the  damping  effect  of 
the  wall  on  the  turbulent  fluctuations.  The  parameter  A  is 
usually  referred  to  as  a  damping  constant.  The  exponential 
term  approaches  zero  at  the  outer  edge  of  the  viscous  sub¬ 
layer  so  that  the  law-of-the-wall ,  equation  (42),  is  valid 
and  it  accounts  for  the  effect  of  kinematic  viscosity  on 
the  turbulence  near  the  wall.  The  damping  constant  A  is 
given  as: 


VTV  V  ^  V  w  V  \~  ^  r  ‘  ;~"  ‘*r  W  r'  *  v 


▼  -'«■ -:  rsr:  r;  v-  r,  i"  r, 


A  =  26v 


ta 


(44) 


with  the  subscript  w  denoting  values  at  the  wall.  Equation 
(43)  was  developed  for  a  fjLat  plate.  To  account  for  flows 
with  non-zero  pressure  gradient,  the  constant  A  is  modified 
as  follows.  From  the  momentum  equation,  it  follows  that 
the  shear  stress  close  to  the  wall  may  be  written  as: 


T  =  T. 


•(Sf> 


(45) 


If  A  were  redefined  to  be  26v(t/p)  "  then  Eq  (45)  leads  to: 


("26vf  +  d£  %  "l-i'l 

L  Lp  dx  p  J  J 


Then,  the  corresponding  expression  for  the  inner-rogion  eddy 


viscosity  becomes: 


xr  2  2 

Z  .  =  K.  y 

inner  1  J 


A  ->  2 


-  exp  -£/  _w  +  dp  z 
.  2ov 1  p  dx  p 


9u I  (47) 
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Outer-Region  Model.  The  eddy  viscosity  in  the  outer 

t 

region  is  given  by  the  expression: 


„/°°U  -u)ay 


outer  “2  o  e 


(0, 


=  K«  pUe  6^ 


(48) 


6  ”  =  ;ye  n  -  u  i dy 


where 


(49) 


In  order  to  account  for  the  intermittent  character  of 


the  outer-layer  flow,  equation  (45)  is  modified  by  an  inter- 
mittency  factor  obtained  by  Klebanoff  (Ref.  28). 


£  =  K2  pUe  <5 

TT  o  u 


(50) 


where  the  transverse  int ermittency  factor  y(y)  is  defined  as: 

y  =  h  1  -  erf{5(y/6  -  0.78)} 
and  app1  -'ximated  by 

Y  =  {1  +  5.5  ( y / 5 ) 6 } ~ 1  (51) 

which  is  a  convenient  and  sufficiently  accurate  approxima¬ 
tion  to  the  error  function. 

The  choice  of  and  Kg  in  the  eddy-viscosity  formulae 
depends  slightly  on  the  definition  of  the  boundary-layer 
thickness  6.  In  several  previous  studies,  (for  example 
Ref.  29:174-191),  the  values  of  the  constants  and  Kg 
are  taken  to  be  0.4  and  0.0168  respectively,  and  5  is 
defined  as  the  normal  distance  from  the  surface  to  a  point 
in  the  field  at  which  F  was  equal  to  0.995. 

The  constraint  used  to  define  the  inner  and  outer 
region  is  the  continuity  of  the  eddy  viscosity.  Starting 
from  the  wall,  the  expression  for  the  inner-eddy  viscosity 
applies  until  =  eq,  provided  that  Y  +  >_50.  If  Y  +  <50  when 
ei  =  eo’  the  switch  point  is  delayed  until  Y+  becomes  50 

and  ei  becomes  equal  to  at  Y+=50.  This  prevents  the 
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suppression  of  the  fully  turbulent  portion  of  the  velocity 
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profile  that  can  occur  at  lower  values  of  R 


This  effect 


"0 


persists  to  higher  and  higher  values  of  R  as  the  mach 

e 

number  increases.  Figure  3  shows  a  typical  eddy-viscosity 
variation  across  the  boundary  layer  for  flow  over  a  flat 
plate . 

Eddy  Conductivity 

The  eddy  conductivity  is  formulated  in  terms  of  a 
static  turbulent  Prandtl  number  P  ,  and  the  eddy  viscosity 

I'  ,  b 

e  {see  equations  (9)  to  (11)}.  The  two-layer  concept  for 

the  eddy  viscosity  model  suggests  that  there  should  also 

be  a  two-layer  model  for  static  turbulent  Prandtl  number. 

Numerous  assumptions  have  been  made  concerning  the  eddy 

conductivity,  and  one  of  the  earliest  assumptions,  which 

has  been  used  extensively,  employed  a  constant  value  of 

unity  for  the  static  turbulent  Prandtl  number.  However, 

experimental  data  definitely  shows  that  P  ,  is  a  function 

r ,  o 

of  (y/6);  hence,  this  assumption  is  expected  to  lead  to 
error . 

The  incompressible  data  indicates  that  P  ,  ranges 

r ,  u 

between  0.7  and  0.9.  Simpson,  Whitten  and  Moffat  (Ref.  30) 
found  that  P  ,  ranges  from  approximately  0.95  at  ( y / 6 ) 

=  0.1  to  0.45  at  (y/d)  =  1.0.  The  data  in  this  region  were 
predicted  well  by  the  expression 


Pr>t  =  0.95(1  -  0 . 5 (y/6 )  } 


as  proposed  by  Rotta  (Ref.  23). 
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For  compressible  flow,  where  very  limited  experimental 
data  are  available,  ^  has  a  value  very  near  unity  in  the 
outer  region  of  the  boundary  layer  and  a  value  between  0.7 
and  0.9  at  the  wall  boundary  (Ref.  24).  Meir  and  Rotta 


(Ref.  25)  found  that  for  1.75  <  M  <  4.5,  P_  .  increased 

00  — •  r  |  y 

above  unity  for  T+<50  and  ranged  between  0.8  and  0.85  as 

the  outer  edge  of  the  boundary  layer  was  approached. 

The  assessment  of  the  current  value  for  P  ,  must  be 

r ,  t 

based  upon  the  agreement  between  experimental  and  calcu¬ 


lated  temperature  profile  over  a  wide  range  of  flow  vari¬ 


ables.  Due  to  inconclusiveness  of  much  of  the  experimental 

data  which  exists  to  date,  a  constant  value  of  P  ,  equal 

r ,  t> 

to  0.9  is  utilized  in  the  present  analysis. 

The  following  chapter,  discusses  the  linearization 
and  discretization  of  the  governing  equation,  followed  by 
the  computational  technique  to  solve  the  system  of  non¬ 
linear  parabolic  finite-difference  equations. 
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Dimensions  and  Distribution 
of  Roughness  Elements  with 
Circular  Cross-section 
( From  Ref  15:8) 


Ill .  Method  of  Solution 


The  system  of  governing  equations  for  compressible 
laminar,  transitional  and  turbulent  boundary  layers  consists 
of  three  coupled  non-linear  partial  differential  equations 
(Equations  (30)  to  (32))  and  two  algebraic  relations 
(Equations  (13)  and  (14)).  The  most  important  feature 
of  this  system  is  that  it  is  parabolic  and,  therefore, 
can  be  numerically  integrated  in  a  step-by-step  procedure 
in  the  streamwise  direction.  In  order  to  cast  the  equa¬ 
tions  into  a  form  in  which  the  marching  procedure  can  be 
efficiently  utilized,  the  derivatives  with  respect  -to  £ 
and  q  are  replaced  by  linear  finite-difference  quotients. 

In  constructing  the  difference  quotients,  the  skenc'-  of 
the  grid-point  distribution  presented  in  Figure  L 
useful  for  reference. 

The  solution  is  obtained  in  the  transformed  plane  for 
arbitrary  grid-point  spacing  in  the  ^-direction  and  for  a 
spacing  in  p-direction  such  that  the  ratio  of  the  spacing 
between  any  two  successive  pairs  of  grid  points  is  a  con¬ 
stant. 

The  advantage  of  having  variable  grid-point  spacing  in 
the  ^-coordinate  becomes  clearly  apparent  for  problems  in 
which  either  the  rate  of  change  of- the  boundary  conditions 
is  large  or  discontinuous  or  the  mean  flow  profiles  are 
changing  rapidly.  Variable  grid-point  spacing  in  the 


31 


o 


^-direction  is  implemented  by  having  very  small  steps  in 
the  initial  region  where  the  flow  gradient  is  very  severe 
and  again  in  some  downstream  region  where  transitional  flow 
exists.  In  the  present  analysis,  certain  stations  in  the 
C-direction  were  designated  where  the  streamwise  step  was 
doubled  from  its  value  at  the  preceding  station.  Hence, 
larger  A £  step  sizes  are  utilized  as  the  flow  progresses 
downstream.  Variable  spacing  in  n -direction  is  even  more 
critical  since  the  flow  gradients  near  the  wall  are  ex¬ 
tremely  large,  whereas  these  gradients  vanish  near  the 
edge  of  the  boundary  layer.  The  relationship  between  At), 
for  the  chosen  grid-point  spacing  is.  given  by  the  following 
equation  (Ref,  19:32-33) 

An  ±  =  (k)i_1  Ar^  (i  =  1 ,  2,  3 - N)  (52) 


where  k  is  the  ratio  of  any  two  successive  steps. 


Difference  Equations 

Three-point  central-difference  relations  in  n-direction 
and  two-point  backward  or  three-point  backward  difference 
relations  in  the  ^-direction  are  used  to  reduce  the  trans¬ 
formed  continuity,  momentum,  and  energy  equations  (Equations 
(30)  to  (32))  to  a  system  of  coupled  difference  equations. 
Following  the  quasilinearization  (see  Appendix  C)  of  the 
non-linear  terms,  the  difference  quotients  produce  linear 
difference  equations  when  substituted  into  the  continuity, 
momentum,  and  energy  equations.  The  resulting  difference 
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difference  equations  are  written  symbolically  as  follows: 
Momentum : 

A1(n,1)Fn_1  +  A1(n,2)Fn  +  A1(n,3)Fn+1 

+  Bl(n,l)en_1  +  B1(n,2)0n  +  B1(nf3)0n+1 
+  CKn.OV^  +  C1(nf2)Vn  +  C1(n,3)Vn+1  =  D1(n)  (53) 

Energy: 

A2(n,1)Fn_1  +  A2(n,2)Fn  +  A2(n,3)Fn+1 

+  B2(n,1)0n_1  +  B2(n,2)9n  +  B2(n,3)9n+1 
+  C2(n,l)Vn_1  +  C2(n,2)Vn  +  C2(n,3)Vn+1  =  D2(n)  (54) 

Continuity: 

+  A3(n,2)Fn  +  A3(n,3)Fn+1 
+  B3(n,1)0n_1  +  B3(n,2)9n  +  B3(n,3)en+1 
+  C3(n,1)Vn_1  +  C3(n,2)Vn  +  C3(n,3)Vn+1  =  B3(n)  (55) 

The  coefficients  A1  1)...,  B1(n,  )...,  Cl(n,1)...,  and 
D1 (n) . . . ,  etc.,  are  functions  of  known  quantities  at 
stations  i-1  and  i-2  and  are  detailed  in  Appendix  D.  The 
dependent  variables  F  ,  and  V  appear  in  a  linear  form 


as  unknown  at  station  i  but  are  assumed  to  be  known  at 
stations  (i-1)  and  (i-2).  Since,  the  equations  comprise 
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a  system  of  parabolic  partial  differential  equations, 
they  can  be  solved  by  inarching  along  £  and  a  Gaussian 
elimination  technique  along  n . 

Computer  Code 

The  computer  code,  ITRACT,  available  in  the  Flight 
Dynamics  Laboratory  computed  the  characteristics  of  lami¬ 
nar,  transitional  and  turbulent  flow  for  either  planar  or 
axisymmetric  flov/  over  smooth  surfaces.  This  code  was 
modified  to  include  the  present  analysis  for  rough  sur¬ 
faces.  To  initiate  this  code,  the  following  quantities 
were  specified  as  inputs: 

y  -  the  ratio  of  specific  heat 

P  -  laminar  Prandtl  number 

r 

turbulent  Prandtl  number 

the  exponent  in  Sutherland’s  viscosity  law 


r,  t 
w 

BO 


ltr 


ratio  of  wall  temperature  (T  )  to  free 
stream  stagnation  temperature  (Tq) 

the  streamwise  location  of  transition  from 
laminar  to  turbulent  flov: 


IDIFF  -  a  flagged  quantity,  tc  specify  the  three- 
point  or  two-point  differencing  scheme  for 
the  streamwise  derivative 

J2DA  -  a  flagged  quantity,  tc  specify  a  2-D  or 
axisymmetric  flov 

The  data  regarding  the  size,  shape  and  density  of  the 
roughness  elements  was  provided  to  the  program  in  sub¬ 
routine  RUFVAR  (Appendix  E),  The  computer  code  provided 
description  of  the  boundary  layer  characteristics.  Some 
of  the  output  of  interest  in  the  computer  code  consists 
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of  the  boundary-layer  profiles  for  Mach  number,  static 
temperature,  velocity,  density,  eddy-viscosity,  enthalpy 
and  pressure.  It  also  provides  values  of  boundary-layer 
thickness,  boundary  layer  momentum  thickness,  the  coeffi¬ 
cient  of  friction,  eddy  viscosity,  and  Stanton  number 
which  is  a  description  of  heat-transfer  at  the  surface. 

The  transformation  from  the  (x,y)  plane  to  the  (f ; , n  ) 
plane  casts  the  boundary  layer  into  a  rectangular  grid  of 
nodes  with  the  surface  of  the  model  located  at  the  level 
j=1,  as  shown  in  Figure  4  where  subscript  (i,j)  refers 
to  £ , n  indices.  The  mesh  spacing  in  the  i,j  direction 
are  AC^  and  Ar)^,  respectively,  which  are  not  constant, 
for  reasons  explained  earlier. 

The  computer  code  solves  the  linear  difference 
equations  (53)  to  (55).  'The  solution  of  this  system  of 
equations  was  determined  by  computing  values  of  F,  6,  and 
V  at  each  of  the  nodes  within  the  grid.  With  all  the 
values  of  these  variables  known  at  station  i-2  and  i-1 , 
the  values  of  F,  0,  and  V  were  solved  at  all  points  j  at 
station  i  using  a  three-point  differencing  scheme  and  a 
Gaussian  elimination  technique.  With  the  boundary-layer 
solution  completed  at  the  station  i,  the  computer  code 
marched  to  station  i+1  in  the  streamwise  direction  and 
the  dependent  variables  F,  0,  and  V  were  determined  in  a 
similar  manner  at  the  preceding  station.  The  entire  pro¬ 
gram  was  therefore  a  sequential  solution  at  a  series  of 
^-stations  from  the  leading  edge  to  the  trailing  edge  of 
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the  model.  Hence,  the  computer  code  followed  the  step- 
by-step  procedure  depicted  in  Figure  4.  With  a  program 
listing  included  in  Appendix  E,  the  main  portion  of  the 
logic  presented  in  Figure  5  required  further  explanation; 
hence.  Appendix  F  was  included  to  discuss  four  important 
portions  of  the  code.  These  included  the  non-dimension- 
alization  of  the  working  variables,  the  computation  of 
eddy  viscosity,  Stanton -number ,  skin-friction  coefficient 
and  the  roughness  variables  f(y),  B(y),  and  A  Fortran 

computer  code  key  is  also  included  as  Appendix  G. 
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IV.  Results  and  Discussions 


The  modified  analysis  and  computer  program  are 
capable  of  predicting  the  effect  of  surface  roughness  on 
turbulent  boundary  layers.  Prior  to  using  it  as  a  valid 
predictive  method,  it  was  necessary  to  establish  the 
authenticity  of  the  original  code.  For  this  purpose, 
results  were  first  obtained  using  the  original  code,  for 
turbulent  flow  over  a  smooth  flat  plate  and  compared  with 
the  experimental  data  of  Dr.  Fiore  at  the  Air  Force 
Wright  Aeronautical  Laboratories  and  Dr.  Donald  Coles, 
at  the  Jet-Propulsion  Laboratory,  California.  Thereafter, 
calculations  were  made  for  the  case  of  turbulent  flow  over 
rough  surface  and  the  computed  boundary-layer  profiles 
for  velocity,  static  temperature,  pitot  pressure,  Mach 
number,  density  and  stagnation  temperature  were  compared 
with  the  corresponding  experimentally  obtained  boundary- 
layer  profiles.  Further,  the  computed  profile  for  the 
rough  surface  were  compared  with  the  computed  profiles 
for  the  smooth  surface,  so  as  to  appreciate  the  effect  of 
surface  roughness  on  the  flow. 

In  these  studies,  a  number  of  Important  assumptions 
were  made.  First,  the  boundary-layer  computation  and  the 
comparison  were  performed  for  flow  over  a  flat  plate; 
hence,  the  effect  due  to  the  stagnation  region  at  the 
leading  edge  and  the  shocking  phenomena  were  neglected. 
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Shapiro  has  alluded  to  the  validity  of  the  assumption  of 
free-stream  conditions  existing  at  some  distance  down¬ 
stream  of  the  leading  edge  of  a  plate  (Fig  2  -21 (c)  and 
the  subsequent  text  (Ref.  3 5 : 1 1 49-1 1 50 ) .  Hence,  free- 
stream  conditions  were  assumed  to  exist  downstream  the 
shock  wave.  Further,  the  angle  of  incidence  of  the 
model  was  assumed  to  be  zero  with  respect  to  the  flow  in 
the  free  stream.  Second,  the  boundary-layer  thickness, 

5,  was  minutely  small  compared  to  the  characteristic 
length  L.  Further,  the  pressure  change  across  this 
boundary  layer  thickness  was  negligible.  Third,  the 
problem  was  limited  to  experimental  cases  where  pressure 
change  along  the  streamwise  direction  was  also  negligible. 

Numerically,  d£  was  considered  zero.  Also,  the  flow  was 
dx 

considered  fully  turbulent  at  the  station  where  the  com¬ 
parison  of  the  results  were  made.  Fourth,  the  flow  was 
considered  inviscid  and  potential  beyond  the  edge  of  the 
boundary  layer.  Finally,  the  Navier-Stokes  equations 
were  simplified  to  the  boundary-layer  equations  to  des¬ 
cribe  the  flow  characteristics  for  y<6  (Ref.  3:117-121). 

In  all  cases  presented  herein,  the  gas  is  taken  to 
be  air  and  is  assumed  to  behave  as  a  perfect  gas  with  a 
constant  ratio  of  specific  heats  (y=  1.4),  a  constant 
Prandtl  number  (P  =  0.72),  and  a  constant  static  turbu¬ 
lent  Prandtl  number  (P  .  =  0.9).  The  molecular  viscosity 

r ,  z 

y  is  evaluated  from  Sutherland’s  viscosity  law.  The 
external  pressure  distribution  used  is  either  experimental 
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or  obtained  from  an  exact  solution  of  the  full  inviseid 
Enter  equations. 


Turbulent  Flow  Over  Smooth  Flat-Plate 

Comparison  With  Dr.  Fiore's  Data  (Ref  37).  For  the 
purpose  of  this  study,  the  results  obtained  from  the  orig¬ 
inal  computer  were  compared  with  Dr.  Fiore's  data  for  tur¬ 
bulent  flow  over  a  smooth  flat  plate.  The  experiment  was 
conducted  at  a  Mach  number  of  5.92  and  a  Reynolds  number 

n 

Rex  of  3.6(10)  .  The  free-stream  stagnation  temperature 

and  stagnation  pressure  were  1098°R  or  1995.38  psia, 

respectively.  The  free-stream  conditions  and  the  data 

of  this  experimental  study  are  tabulated  in  Table  I. 

The  computer  code  results  were  obtained  for  the  same 

free-stream  conditions.  The  point  at  which  the  transition 

from  laminar  flow  to  turbulent  was  initiated,  was  obtained 

6 

from  the  transition  Reynolds  number  of  1.5(10)  .  The  wall 
temperature  was  602°R.  The  length  of  the  model  was  17.1 
inches,  and  the  boundary  layer  profiles  computed  at  the 
end  of  the  flat  plate  (i.e.,  X/L  =  1.0)  were  compared 
with  the  experimental  profiles. 

There  is  excellent  agreement  between  the  computed 
and  the  experimental  boundary-layer  profiles  as  shown  in 
Figures  7  through  12,  except  in  the  region  close  to  the 
wall.  The  velocity  profiles  as  shown  in  Fig.  7  compared 
very  well,  except  in  the  sub-layer  region  where  the 
computer  code  overpredicts  the  velocity  slightly.  On 
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the  other  hand,  the  temperature  profile  (Fig.  8)  is  slightly 
underpredicted  in  the  same  region,  which  is  quite  under¬ 
standable  due  to  the  fact  that  the  stagnation  temperature 
remains  constant.  The  same  trend  is  observed  in  the 
boundary-layer  profiles  for  density,  Mach  number,  total 
temperature  and  pressure  is  depicted  in  Figures  9,  10, 

11,  and  12,  respectively.  The  discrepancy  in  the  sub¬ 
layer  region  could  be  associated  with  either  of  the  fol¬ 
lowing  two  sources  of  error: 

(a)  The  Cebeci-Smith-Mosinki 1 s  model  used  in  this 
study  to  describe  the  sub-layer  region  of  the  boundary 


layer  was  not  adequately  taking  into  account  the  inter¬ 
mittent  nature  of  the  turbulence  in  this  region. 

(b)  The  boundary-layer  probe  installed  at  the  base 
of  the  model  could  influence  the  force -balance  measure¬ 
ments.  This  could  result  in  a  lower  experimental  pressure 


ratio  near  the  wall,  thereby  causing  some  error  in  the 
temperature  profiles  (Ref  39). 

The  tvo  possibilities  mentioned  above  were  examined 
further.  It  was  observed  that  the  Cebeci-Smith-Mosinski ' s 
model  has  been  used  extensively  in  many  previous  studies 
(Ref.  19)  and  has  been  found  to  predict  sub-layer  region 
very  accurately.  While  investigating  the  classical  eddy 
viscosity  models  in  an  attempt  to  establish  a  benchmark 
for  future  development  of  turbulent  boundary  layer  re¬ 
search,  Shang,  Hankey,  and  Dwoyer  (Ref.  19)  found  that 
the  Cebeci-Smith-Mosinski ' s  model  was  adequate  to  predict 
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the  boundary  layer  in  the  sub-layer  region.  Hence,  the 
first  possible  source  of  error  was  discarded.  The  size 
of  the  probe  used  in  the  experimental  study  was  of  the 
order  of  the  sub-layer  thickness  which  could  influence 
the  flow,  thereby  leading  to  a  lower  experimental  pitot 
pressure  ratio  near  the  wall  (Ref.  39). 

To  correctly  establish  the  cause  of  the  error  in  the 
sub-layer  region,  it  was  decided  to  carry  out  another 
comparison  of  the  computed  velocity  and  Mach  number 
profiles  with  the  experimental  data  of  Dr.  Cole  (Ref.  38) 
for  turbulent  flow  over  a  smooth  flat  plate. 

Comparison  With  Cole’s  Data  (Ref.  38).  The  comparison 
with  Cole's  data  was  carried  out  to  check  the  accuracy  of 
the  original  code,  since  the  comparison  with  Dr.  Fiore's 
data  showed  some  discrepancy  between  the  computed  and 
the  measured  results  in  the  sub-layer  region.  The  data 
of  Dr.  Cole,  as  tabulated  in  Table  II,  was  obtained  from 
experimental  investigation  of  turbulent  flow  over  a 
smooth  flat-p]ate  for  the  test  conditions  as  follows: 

=  4.554 

P0  =  8132.788  I6f/ft2 
Tq  =  522.64°R 

VTo  =  *676^ 


The  length  of  the  flat-plate  was  21.48  inches  and 
the  boundary-layer  profiles  were  compared  at  X/L  =  1.0. 
In  Figs  13  and  14*  the  comparison  between  the  computed 
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and  experimental  velocity  and  Mach  number  profiles  is 
presented.  It  is  observed  that  excellent  agreement 
exists  between  the  computed  and  experimental  profiles 
everywhere,  including  the  sub-layer  region  near  the  wall. 

Via  this  comparison,  the  authenticity  of  the  origi¬ 
nal  computer  code  was  established.  The  next  step  was  to 
verify  the  roughness  model  and  the  modified  computer  code 
by  conducting  a  comparison  of  the  boundary-layer  profiles 
predicted  by  the  modified  computer  code  with  the  experi¬ 
mental  data  of  Dr.  Fiore  (Ref.  20)  for  flow  over  a  rough 
flat  plate. 


Turbulent  Flow  Over  Rough  Flat-Plate 

Comparison  With  Dr.  Fiore’s  Data  (Ref.  40)  After  the 
authenticity  of  the  original  code  was  established,  the 
next  step  was  to  modify  the  computer  code  and  ther.  verify 
its  results.  Prior  to  comparing  the  results  of  the  modified 
computer  code  with  the  experimental  data  for  turbulent 
flow  over  a  rough  flat-plate,  results  for  turbulent  flow 
over  a  smooth  flat-plate  were  obtained  from  the  modified 
code,  by  specifying  roughness  elements  in  the  code  to  be 
of  zero  dimension.  The  results  obtained  were  in  complete 
agreement  with  those  from  the  original  code  discussed  in 
the  previous  section  and  plotted  in  Figs  7  through  12. 

This  showed  that  the  inclusion  of  the  roughness  model  to 
the  code  has  not  affected  the  computer  code  for  turbulent 
flow  over  smooth  surfaces. 
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Next,  to  verify  the  roughness  model  used  in  this  study, 
the  results  computed  by  the  modified  code  were  compared  with 
the  experimental  data  of  Dr.  Fiore  (Ref  40)  of  AFWAL,  for 
turbulent  flow  over  a  rough  flat  plate. 

The  experiment  was  conducted  at  the  following  test 
condit.  .  - 


M0  =  5.534 

Po  =  2003.53  psia 

T  =  1 1 22. 94°R 
o 

T  /T  =  .5829 
o  w 

Re„  =  3.665(1 0) 7 


The  roughness  element  distributed  over  the  flat  plate 
had  the  following  dimensions  (Fig.  6) 


(i)  Height  of  the  element  =  k  =  0.02  inches 

(ii)  Breadth  of  the  element  =  b  =  0.04  inches 

(iii)  Depth  of  the  element  in  the  flow  direction 
-  0.04  inches 


The  density  of  the  roughness  elements  was  specified 
as  follows  (Fig.  6)  ’• 

(i)  Spacing  between  two  adjacent  elements  loca  sd 
along  a  line  perpendicular  to  the  flow  direction  =  B  = 
.08  inches. 

(ii)  Spacing  between  two  adjacent  elements  located 
along  a  line  parallel  to  the  flow  direction  =  C  =  .08 
inches . 


V  V 
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The  length  of  the  flat  plate  was  17.15  inches,  and 
the  boundary-layer  profiles  computed  at  X/L  =  1,0  were 
compared  with  the  experimentally  obtained  boundary-layer 
profiles.  The  free-stream  conditions  and  the  experimental 
data  are  tabulated  in  Table  III. 

The  results  presented  here  are  quite  encouraging, 
in  that  a  rather  basic  model  yields  results  that  are  in 
good  areeement  with  many  of  the  observed  trends  regarding 
the  influence  of  surface  roughness.  Figure  15  shows  the 
excellent  agreement  obtained  for  the  velocity  profiles. 

The  computed  profile  shows  the  expected  increase  in  thick¬ 
ness  and  change  of  shape  due  to  the  presence  of  the 
roughness  elements.  The  maximum  error  in  the  computed 
and  the  experimental  velocity  is  about  1.88  percent  and 
occurs  near  the  wall,  and  is  attributed  to  the  same  cause 
as  for  the  smooth-wall  case  as  discussed  earlier  in  this 
chapter. 

The  predicted  and  the  measured  temperature  profiles 
shown  in  Fig.  16  show  the  same  trend,  except  that  the 
computer  code  has  overpredicted  the  static  temperature. 

The  distributed  roughness  elements  over  the  surface 
cause  the  flow  to  slow  down  due  to  form  drag  in  the  wall 
region  of  the  elements.  This  results  in  a  higher  static 
temperature  in  order  to  keep  the  total  temperature  the 
same . 

The  initial  rise  in  the  static  temperature  to  about 
827°R  occurs  in  the  region  corresponding  to  the  height  of 


the  roughness  element,  i.e.,  for  y  <_  .02  inches.  At  y  >  h 
(height  of  the  element),  the  flow  is  no  longer  obstructed 
by  the  roughness  elements  and  the  temperature  profile 
follows  the  same  decreasing  trend  as  that  for  the  smooth 
surface  profile.  In  the  experimental  results,  the  static 
temperature  is  approximately  670°R  at  y  =  .02  inches  and 
the  error  is  about  20$  l.hlu  point.  This  discrepancy 
occurs  in  the  sub-layer  region  and  could  be  traced  back 
to  the  influence  of  the  probe  on  the  flow  (Ref  39).  The 
thin  wall  region  plays  a  very  dominant  role  in  determining 
the  entire  structure  of  the  boundary-layer. 

The  density  profiles  shown  in  Fig.  17  depicts  a 
trend  consistent  with  that  of  the  static  temperature 
profiles  since  the  density  at  a  point  in  the  boundary- 
layer  is  simply  the  reciprocal  of  the  corresponding 
static  temperature.  The  profiles  are  in  reasonably  good 
agreement,  except  in  the  region  close  to  the.  edge  of  the 
boundary  layer. 

The  experimental  boundary-layer  thickness  (6)  is 
.383  inches  and  the  measured  data  exhibits  a  r..onotonic 
increase  in  the  density  from  the  wall  to  the  edge  of  the 
boundary  layer.  The  same  trend  is  observed  in  the  com¬ 
puted  profile  also.  Hence,  if  the  density  were  plotted 
against  y/6,  for  0  <  y/6  <  1,  then  the  profiles  would 
appear  to  be  in  reasonable  agreement.  As  y  increases 
slightly  beyond  0.4  inches,  there  is  a  sudden  decrease 
in  the  experimental  density  values,  with  the  free-stream 
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density  value  being  attained  shortly  beyond  y  =  0.4. 
Actually  this  trend  is  not  consistent  with  the  static 
temperature  profile  given  in  the  experimental  data 
(Fig.  16).  . 

Comparison  of  Computed  Rough  and  Smooth  Boundary  Layer 
Results 

The  velocity  profiles  with  and  without  roughness  are 
compared  in  Fig  21 .  Both  the  profiles  are  for  turbulent 
flow  over  a  flat  plate.  The  two  cases  have  slightly 
different  free-stream  conditions,  as  tabulated  in  Table 
IV.  The  smooth-wall  profile  exhibits  fully  developed 
turbulent  flow  characteristics  and  a  smooth  monotonic 
increase  in  velocity  from  zero  at  the  wall  to  the  free- 
stream  value  at  the  edge  of  the  boundary  layer. 

The  effect  of  distributed  surface  roughness  on  the 
boundary  layer  development  can  be  ovserved  as  a  reduction 
in  the  values  of  the  velocity  in  the  wake  of  the  roughness 
elements.  In  the  roughness  model,  the  rough  surface  is 
idealized  as  being  made  up  of  identical  elements.  The 
form-drag  description  of  the  roughness  element  causes  the 
flow  to  slow  down  in  their  wake  region;  hence,  the  rough- 
wall  profile  is  not  as  full  as  that  corresponding  to  the 
smooth  wall.  Table  V  shows  a  corresponding  increase  in 
all  the  boundary-layer  integral  thicknesses  when  roughness 
is  included. 

Semi-empirical  analysts  of  turbulent  boundary  layers 


are  often  based  on  some  form  of  law-of -the-wall  correlation, 
The  law-of-the-wall  is  defined  most  simply  in  terms  of  the 


dimensionless  parameters 


U*  =  U  and  Y+  =  YIIt 


where 


/T  /p  ^ 
\  w'  Kw  / 


(56) 


(57) 


Within  the  boundary  layer,  three  distinct  regions  are  found 
to  exist;  the  sub-layer,  the  logrithmic  region  and  the 
velocity-defect  region.  The  general  form  of  the  relation¬ 
ships  governing  smooth-wall  boundary  layer  is  given  by: 


U+  =  Y+ 


U  |  In  I*  t  C 


for  Y  <  11 


for  Y  >  11 


(58) 


(59) 


The  effect  of  roughness  on  the  law-of-the-wall  equation 
has  been  shown  uo  result  solely  in  a  shift  in  the  inter¬ 
cept,  C,  in  equation  (59).  The  same  trend  has  been 
predicted  in  the  computed  smooth  and  rough  velocity 
profiles.  The  rough-wall  profile  looks  similar  to  the 
corresponding  smooth-wall  profile,  except  that  it  is 
shifted  towards  lower  values  of  velocities.  This  trend 
is  consistent  with  the  rough-wall  results  where  the  shift 
in  the  profile  is  directly  related  to  the  size  and  density 
of  the  roughness  elements.  The  computed  boundary-layer 
profiles  show  the  expected  increase  in  thickness  and 
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change  of  shape  due  to  roughness 


Skin  Friction  and  Heat-Transfer 

The  addition  of  roughness  elements  to  the  smooth 
surface  can  significantly  affect  turbulent  skin  friction 
and  heat  transfer.  The  roughness  model  used  in  this 
study  was  tested  for  prediction  of  skin  friction  augmen¬ 
tation  due  to  addition  of  roughness  element  to  a  smooth 
surface.  A  comparison  between  the  computed  skin  friction 
values  for  a  smooth  surface  and  a  rough  surface  was 
accomplished.  The  predicted  skin  friction  coefficient 
for  the  two  cases  was  plotted  versus  dimensional  distance 
along  the  flat  plate  in  Fig  29.  In  this  study,  square 
roughness  elements  with  K  =  .02"  b  =  d  =  .04"  and  B  =  C  = 
.08"  were  used.  Fig  29  shows  the  expected  increase  in 
the  skin  friction.  For  rough  surface  the  skin-friction 
was  40-60$  higher  as  compared  to  smooth  surface.  The 
validity  of  the  skin  friction  augmentation  by  the  model 
could  not  be  confirmed  due  to  non-availability  of  any 
experimental  data. 

As  regards  heat  transfer,  the  original  model  for 
calculation  of  Stanton  number  was  not  modified  in  this 
study,  hence,  the  expected  increase  in  the  heat  transfer 
due  to  roughness  elements  was  not  predicted  by  the  code. 
On  the  other  hand,  a  decrease  in  Stanton  number  was  noted 
since  smooth  surface  had  a  higher  pressure  gradient 
compared  to  the  rough  surface  near  the  wall.  Hence,  a 
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need  for  incorporation  of  a  suitable  model  for  accurate 
prediction  of  heat  transfer  augmentation  was  felt.  It  may 
oe  noted,  that,  one  particular  useful  aspect  of  this  rough 
wall  turbulence  model  is  that  the  results  can  be  examined 
to  determine  the  nature  of  the  roughness  influence  on 
turbulent  boundary-layers.  One  rather  conspicuous  con¬ 
clusion  is  that  the  Reynolds  analogy  between  friction  and 
heat  transfer  is  not  preserved  with  significant  roughness. 
This  result  is  well  known  and  derives  from  -.he  absence  of 
a  heat  transfer  analogy  to  form  drag  on  elements.  The 
computation  shows  that  the  velocity  fluctuations  increase 
in  proportion  to  friction  velocity  lT  +  =  (Tv/o  )i'  but  the 
temperature  fluctuations  are  hardly  changed  by  roughness. 
Since  T  a-  u  /  v  '  and  q  a  v’  tT  * ,  the  heat  transfer  augmenta- 
tion  is  the  square  root  of  the  skin  friction  augmentation: 

St  1  =  Gf  <  f 

St  "  1  Cf  “  1  (60) 

0  0 

where  subscript  Tof  denotes  smooth  wall. 
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Fig,  8,  Comparison  of  Computed  Temperat' 
Profiles  for  Smooth  Surface  Wit 
Experimental  Data  of  AFWAL 
(Ref  37) 
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Fig.  10.  Comparison  of  Computed  Mach 
Number  Profile  for  Flow  Over 
Smooth  Surface  With  AFWAL 
(Ref  37) 
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15.  Comparison  of  Computer  Velocity 
Profile  for  Rough  Surface  With 
Experimental  Data  of  AFWAL) 

(Bef  40) 
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7.  Comparison  of  Computed  Density 
Profile  for  Rough  Surface  With 
AFWAL  (Ref  40) 
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Fig.  18,  Comparison  of  Computed  Mach 
Number  Profiles  for  Rough 
Surface  With  AFWAL  (Ref  40) 


44  4 


Pressure  (psi) 


Fig.  20.  Comparison  of  Computed  Pressure 
Profile  for  Rough  Surface  With 
AFWAL  Data  (Ref  40) 
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Fig.  21 .  Comparison  of  Computed  Velocity 
Profiles  for  Flows  Over  Smooth 
and  Rough  Surfaces 
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Table  I 


Data  for  Smooth  Flat  Plate  (Ref  37)  Probe 
Boundary  Layer  Profiles 


Y 

Inch 

Pt(y) 

psia 

Tt(y)  + 
°R 

M(y) 

T  (y ) 

°R 

U(y) 
ft/ sec 

f  (y)xio4 
jj  sec^/ft 

0 

1  .385 

604.62 

0 

604.62 

0 

1  .9234 

.021 

5.65 

963.35 

1 . 666 

619.37 

2032,11 

1 .8776 

.053 

12.52 

1034.0 

2.578 

443.87 

2661 .67 

2.6199 

.083 

19.98 

1066.0 

3.292 

336.54 

2959.25 

3.4555 

.087 

20.73 

1069.0 

■  3.355 

328.75 

2981 .07 

3.5374 

.112 

26.25 

1089.0 

3.789 

281  .28 

3113.96 

4-1345 

.138 

32.01 

1094-0 

4.194 

242.12 

3197.95 

4. 8031 

.148 

35.75 

1 097.0 

4.438 

222.13 

3240. 80 

5.2352 

.169 

42.39 

1099.0 

4.840 

193.35 

3297.33 

6.0147 

.204 

52.14 

11 00.0 

5.375 

162.27 

3355.21 

7.1667 

.205 

53.94 

1100.0 

5.469 

157.57 

3363.61 

7.3804 

.232 

60.53 

1100.0 

5.797 

142.47 

3390.^6 

8.1 627 

.255 

61 .23 

1100.0 

5.831 

141 .03 

3393.00 

8.2458 

.272 

63.50 

1100.0 

5.939 

136.57 

3400.88 

8.5153 

.284 

62.73 

1100.0 

5.903 

138,05 

3398.27 

8.4239 

.322 

63.18 

1100.0 

5.924 

137.18 

3399.80 

8.4773 

.328 

63.00 

1100.0 

5.915 

137.53 

3399.19 

8.4560 

.359 

63.10 

1100.0 

5.920 

137.33 

3399.53 

8.4678 

.378 

63.31 

1100.0 

5.930 

136.93 

3400.24 

8.4928 

.395 

63.09 

1100.0 

5.920 

137.35 

3399.50 

8.4666 

.411 

62.90 

1100.0 

5.911 

137.72 

3398.85 

8.4441 

.445 

63.37 

1100.0 

5.933 

1 36.82 

3400.45 

8.4999 

.447 

62.97 

1100.0 

5.914 

137.58 

3399.09 

8.4524 

.476 

63.18 

1100.0 

5.924 

137.18 

3399.80 

8.4773 

.503 

63.91 

1100.0 

5.959 

135.79 

3402.25 

8.5640 

.512 

63.18 

1100.0 

5.924 

137.18 

3399.80 

8.4773 

.535 

63.50 

1100.0 

5.939 

136.57 

3400.88 

8.5153 

.565 

63.58 

1100.0 

5.943 

136.42 

3401  .1  5 

8.5248 

.569 

63.65 

1100.0 

5.946 

136.28 

3401 .39 

8.5331 

.596 

63.17 

1100.0 

5.924 

137.20 

3399.77 

8.4761 

.631 

63.24 

1100.0 

5.927 

137.06 

3400.01 

8.4845 

Table  II 


Data  of  Dr.  Cole's  Experiment  for  Turbulent 
Flow  Ove-^  Smooth  Surface  (Ref  38) 


111 

u/u^ 

U £ 

M/M 

0 

0.60 

0.500 

0 . 60 

,275 

0.80 

0.555 

1.19 

.355 

0.98 

0.60 

1.60 

.385 

1.50 

0.645 

2.20 

.412 

2.10 

0.687 

2.80 

.439 

3.80 

0,750 

3.70 

.480 

5.20 

0.790 

5.50 

.535 

6.10 

0.81  5 

6.65 

.575 

6.80 

0.030 

7.85 

.61  0 

7.85 

0.845 

9.45 

.659 

9.45 

0.880 

10.70 

.700 

10.55 

0.900 

12.30 

.756 

13.60 

0.950 

13.65 

.  799 

15.60 

0.975 

15.70 

.  880 

16.25 

0.985 

17.19 

.939 

19.20 

1  .00 

19.30 

.980 

21.20 

1.00 

21  .20 

.992 

22.80 

1  .00 

22.80 

1.00 

24.75 

1.00 

24.75 

1  .00 
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Table  III 

Data  for  Rough  Flat  Plate  (Ref  40)  Probe 
Boundary  Layer  Profiles 


rt(y) 

psia 


Tt(y) 


M(y) 


T(y) 

°R 


u(y) 

ft/sec 


f (v)x104 
#sec^ / ft4 


0 

1.632 

654.52 

0 

654.52 

0 

2.0936 

.024 

4.97 

925.0 

1.3992 

664.72 

1767.63 

2.061  4 

.065 

9.19 

1 000.50 

1.9980 

558.84 

2314.34 

2.4521 

.092 

10.93 

1 035.0 

2.1 962 

526.80 

2470.01 

2.6012 

.103 

13.05 

1045.0 

2.41 56 

482.22 

2599.26 

2.8417 

.126 

14.34 

1 060.0 

2.5398 

462.86 

2677.44 

2.9605 

.147 

16.50 

1 072.0 

2.7351 

429.47 

2777.34 

3.1907 

.171 

20.61 

1085.0' 

3.0724 

375.70 

2918.07 

3.6473 

.181 

20.50 

1089.0 

3.0639 

378.46 

2920.62 

3.6208 

.212 

26.28 

1 098.0 

3.4845 

320.27 

3055.59 

4.2786 

.230 

27.50 

1104.0 

3.5670 

311.46 

3084.56 

4.3997 

.246 

30.05 

1106.0 

3.7334 

292.00 

3126.03 

4.6928 

.267 

34.80 

1110,0 

4.0251 

261 .77 

3191 .08 

5.2348 

.297 

44.20 

1114.0 

4.5476 

216.89 

3281 .73 

6.31 79 

.299 

44.10 

1115.0 

4.5424 

217.49 

3282.47 

6.3005 

.332 

55.70 

1116.0 

5.1147 

1 79.07 

3353.78 

7.6522 

.355 

62.50 

1117.0 

5.4222 

162.35 

33.85.34 

8.4404 

.364 

64.09 

1117.0 

5.4917 

158.85 

3391 .54 

8.6263 

.383 

65.08 

1117.0 

5.5344 

1 56.75 

3395.26 

8.7420 

.412 

67.25 

1117.50 

5.6271 

1 52.40 

3403.83 

8.9917 

.426 

67.60 

1117.50 

5.6419 

1  51 .71 

3405.05 

9.0326 

.447 

67.50 

1117.50 

5.6377 

1  51 .90 

3404.70 

8.6703 

.475 

64.50 

1117,50 

5.5094 

1 58.05 

3393.85 

8.6002 

.480 

63.90 

1117.50 

5.4834 

1  59.34 

3391 .57 

8.4819 

.503 

62.92 

1118.0 

5.4407 

161.56 

3388.52 

8.5601 

.530 

63.59 

1118.0 

5.4699 

160.08 

3391  .1  4 

8.6704 

.547 

64.60 

1119.0 

5.5137 

1 58.05 

3396.51 

3.5968 

.565 

63.97 

1119.0 

5.4864 

1 59.40 

3394.12 

8.5968 

.590 

63.77 

1119.0 

5.4773 

159.83 

3393.35 

8.5735 

.607 

64.50 

1119.0 

5.  5094 

1 58.26 

3396.13 

8.6587 

.670 

64.27 

1119.0 

5,4995 

158.75 

3395.26 

8.6319 
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Table  IV 

Free  Stream  Condition  for  Smooth  and  Rough 
Experimental  Data  (Ref.  37>  40) 


Smooth  Surface 


5.92 


Rough  Surface 


5.534 


To  (0R) 


Po  (psia) 


Tw/To 


L  (in) 


1087 


1980 


1122.94 


2003.53 


.  5538 


.  5829 


17.15 


17.20 


Table  V 


Comparison  of  Thickness  for  Smooth  and  Rough  Cases 


Parameter 

Smooth 

Rough 

Me 

5.92 

5.54 

P0(psia) 

1995.38 

2003.53 

W 

1098.38 

1122.94 

.5538 

.5884 

CO 

CD 

16701 .92 

62954.7 

6  (in) 

.2069 

.  4735 

6  * ( in ) 

.0992 

.3792 

6  (in) 

.00793 

.0295 
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V.  Conclusions  and  Recommendations 


The  original  computer,  ITRACT,  solved  the  compres¬ 
sible  turbulent  boundary  layer  over  smooth  surface.  The 
purpose  of  this  study  was  to  extend  the  usefulness  of  the 
computer  code  by  incorporating  a  modification  for  inclu¬ 
sion  of  the.  effect  of  surface  roughness  on  compressible 
turbulent  boundary  layer.  In  this  study,  the  roughness 
model  proposed  by  Finson  and  Clark  and  followed  by 
Christoph  and  Fletcher  was  employed  without  invoking 
the  modification  of  the  turbulence  model.  Roughness 
was  represented  by  distributed  sources  and  sinks  in  the 
appropriate  governing  equations.  The  most  important  term 
was  a  sink  term  in  the  mean  momentum  equation  representing 
form  drag  due  to  roughness  elements.  The  governing 
boundary  equations  were  cast  in  a  form  to  account  for 
blockage  effects  of  the  roughness  elements.  With  the 
computer  code  modified  to  include  the  effect  of  surface 
roughness,  compressible  turbulent  boundary  layer  flows 
perturbed  by  the  roughness  elements  could  be  solved. 

The  roughness  model  employed  and  the  modified  computer 
code  were  verified  through  comparison  of  the  calculated 
results  with  the  experimental  data  of  Dr.  Fiore  for 
compressible  turb  llent  flow  over  rough  surface. 

As  a  result  of  incorporating  a  roughness  model  in 
the  computer  code  to  predict  the  effect  of  surface 
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roughness  on  the  compressible  turbulent  boundary-layer,  the 
following  conclusions  were  reached: 

1.  The  results  presented  here  are  quite  encouraging, 
in  that  a  rather  basic  model  yields  results  that  are  in 
agreement  with  many  observed  trends  regarding  the  influence 
of  surface  roughness. 

2.  The  assumptions  inherent  to  this  model  are  limited 
to  the  basic  nature  of  the  flow  around  the  roughness  ele¬ 
ments,  and  r.c  approximations  have  been  made  regarding 
profiles  of  the  boundary  layer  quantities,  turbulence 
level,  or  relations  between  the  momentum  and  energy 
fluxes . 

3.  The  computer  code  is  capable  to  handle  roughness 
elements  cf  any  general  shape.  It  is  not  restricted  to 
roughness  elements  with  circular  cross-section  only. 

The  following  areas  are  identified  for  further  study  and 
improvement  to  the  method  described  and  the  -computer  code: 

1.  The  Levy-Lees  transformed  coordinate  system  could 
not  properly  capture  the  turbulent  boundary  layer  thick¬ 
ness,  specifically,  for  the  rough  surface  case.  The  trans¬ 
formed  cocrd.inaie  q  became  too  large  and  it  was  necessary 
to  monitor  the  numerical  solutions  and  add  points  in  the 
outer  region  tc  accommodate  the  boundary  layer  growth. 

A  better  turbulence  grid  generation  for  this  type  of  flow 
can  be  used,  such  as  that  proposed  by  Carter  et  al. 

2.  The  form  drag  coefficient  and  spacing  between  the 
elements  be  allowed  to  vary. 


v  .  v*  rvr:  r: 


i 

3.  The  effect  of  one  element  on  another  be  included.  ‘ 

I 

I, 

4.  It  would  be  useful  to  have  better  data  to  study  s 

i 

on  the  effects  of  roughness  density  and  shape  for  distri¬ 
buted  roughness.  •' 

» 

5.  The  combined  effect  of  roughness  and  mass  addition  J 

and  the  effect  of  roughness  at  strongly  supersonic  or  \ 

hypersonic  edge  Mach  numbers  may  be  examined. 

«■ 

6.  A  model  to  accurately  predict  the  heat  transfer  3 

augmentation  due  to  roughness  elements  be  incorporated.  *! 

kj 

Moreover,  the  calculated  skin  friction  augmentation  due  ‘j 

to  roughness  elements  be  compared  with  experimental  data  | 

to  establish  its  accuracy.  !j 

■1 

d 

3 
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Appendix  ’A1 

Derivation  cf  Source  Term  for  Energy  Equation 


Including  a  sink  term  R  =  -1_  pu  un  D(y )  in  the  mean 

u  2  u  2 

momentum  equation  (equ.  )  requires  inclusion  of  an  appro- 

Q 

priate  source  term  =  1_  pu  D(y )  in  the  static  enthalpy 

^  a2 

equation  such  that  the  total  enthalpy  is  not  altered. 

The  mean  static  enthalpy  equation  after  addition  of 


source  term,  R^,  may  be  written  as 

pu  3_h  +  pu  3h_  =  u  djq  +  3_  j  u_  3_ 
3x  3y  dx  3y  ^  Pr  3 


3h  -  ph ' v  1 

3y 


(U  3jJ  -  puv  1  3_u 

3y  /  3y 


+  1  Cp.  pu3  D(y ) 

r\ 


(K1) 


Rewriting  the  mean  momentum  equ 


dp  =  3 
dx  3y 


U  _3u 

ay 


p  u'v 


-  p  u  3  u  4  p  v  3  u 
3x‘  3y 


1  pu  * 

2 


2lii 

£2 


(A. 2) 


Substitute  (A. 2)  in  (A.1)  to  get, 

P  u/  3_h  +  u  3_u\  +  pv  /  3h_  +  u  3_u  | 
y  3x  dx  j  y 3y  3y  / 


+  u  3ju  I  +  pv 
3x  / 


r  \  *  u  f 

d_  j  u  u  3u_  -  u  pu'TT-  +  u_  _3h  -  ph '  ~  ( A .  3 

3y  !  3y  Pr  3y 

in  term  of  total-enthalpy,  the  fluctu- 


To  express  ph'v',  in  term  of  total-enthalpy,  the  fluctu- 

ay 

ating  total  enthalpy  may  be  expressed  as:  (Ref  3c;76~77) 

H«  =  h’  +  u  u>  (A. 4) 
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Multiplying  both  sides  (A. 4)  by  pv*  and  averaging, 


p  H'v’  =  p  h '  v '  +  u  p  u'v 


(A. 5) 


Noting  that, 


3H  =  3h  +  u  3u 
8x  3x  3x 


(A. 6) 


and. 


3H  =  3h  +  u  3u 
3y  3y  3y 


(A. 7) 


Substituting  equations  (A. 5),  (A. o',  and  (A, 7)  in  (A. 3) 


pu  aH  +  pv  aH  =  a_ 
3x  3y  3y 


u  u  a_u  -  c '  v  > 

3y 

+  u_  3_H  -  ^  u  3u 
Pr  3y  Pr  ay 


(A. 8) 


Heat  transfer  normal  to  the  mair.  11;-.:,  |y,  is  given  by 
(Ref  36:66) 


-v 


n_  si 

Pr  3y 


(A. 9) 


Substitute  (A. 9)  in  (A, 8)  to  get, 

pu  aH  +  pv  3H  =  a_  up  -  9y  -  pH'v1 
3x  3y  3y 


(A. 10) 


It  is  evident  from  equation  (A.10),  which  is  the  total- 
enthalpy  equation,  that  the  addition  of  the  source  term,  R^, 
to  the  static  enthalpy  equation,  has  not  altered  the  total 
enthalpy  equation. 


I 


Appendix  1 B  * 

Derivation  of  Conservation  Equations  Including 
the  Roughness  Effect  in  Levy-Lees  Variables 


(B1 )  Roughness  Effects 

Equations  (15),  (16),  (20),  and  (21)  are  repeated  here 
for  convenience: 


Sink  term  included  in  the  mean  momentum  equation: 

Ru  =  -1  pu2  C_  Dili 


Source  term  included  in  the  mean  static  enthalpy  equation: 

Rh  --  +1  pu3  CD  D(xl 

£ 


( B .  1  ) 


(B.2) 


Blockage  terms: 


B(y )  = 


f(y)  = 


^1  -  /  B  ( y ) 


(B.  3 ) 


(B.4) 


( 32 )  Conservation  Equations  in  Physical  Plane 
(a)  Continuity  Equation: 

Introducing  the  blockage  effect  in  continuity  equation 

as  explained  earlier  in  Chapter  II,  Equation  (4)  becomes 

1  -  D(y )  3_  (r^pu)  +  1  -  TiP^(y )  9  (r^'pv)  =  0  (B.4) 

£  3x  / n  2  3y 


Dividing  thru  by  B(y),  (B4)  becomes: 

f(y)  3_  (r^pu)  +  3 _  (r^pv)  =  0 

3x  3y 


( B.  5 ) 


A  - 


v-‘* 


■MITPh'V’ 

■;3 


r  r.  Mv 


V.v.v/ 


vv.vs1 


V  \ 

ft** ' 

AV 


(b)  Momentum  Equation: 

Introducing  the  blockage  effect  and  adding  the  sink  term, 
the  momentum  equation  (equation  (5))  becomes: 
pu  j  1  -  D(y)  (  3u  +  pvj  1  -  ttD2(y)  I  3u  - 

j  t  Jrc  j  — £r[s7 


. .  im(  u .  r .  i_ 


3  r^  (p3_u  -  pu  ’  v  *  ) 
3y  3y 


-  1  pu  cn  D(y)  *(B.6) 

r~s  XJ  r\ 


O 


Dividing  thru  by  B(y),  ( B , 6 )  becomes: 

f(y)  pu  3u  t  pv  3u  = 

3x  3y 

-f(y)  d£  +  1_  1 _  3_ 

dx  J  B(y)  3  y 


B(y  )r^  ( p3u_  -  pu  1  v  '  ) 

3y 


D(y) 


-  1  PU2  CD_ 

2  B(y)  r 


( B.  7 ) 


( o )  Energy  Equation : 

Introducing  the  blockage  effect  and  adding  the  source 
term,  the  energy  equation  (equation  (6))  becomes: 


D(y)  (  pu  3 ( CpT  )  +  F ( y )  ov  3_  (CpT) 


3x 


9y 


u 


+  B(y)  u 


R(zl  l  d£  +  1_  jL.  B(y)  1'^  iii  3  (CpT)  1 
*  (  dx  rj  3y  |  Cp  3y  j 

/  3u  \  2  +  L  Li"  B(y)rj  (-Cpp7rTT)l 

V 3y  I  P  3y  L  J 

-  B  ( y )  ou  '  v  '  3_u  +  1_  pu^  Cn  D(y 


2  L2 


(B.8) 


<> 

H 
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L  IjL  1*1  f  1  M  I  tt.  &  >  »'►  f  -t  mk x\j  k'  ,  *»*■ 


|»*»\«%  V' 


rv*T  v^’v.'nvvi  v ^ I'v  ■%-  v  vi  v.  svw 


„-  *,-  v  '.'.vviv,  yy«t 


0 


Dividing  thru  by  B(y),  Equation  (B.6)  becomes: 

f (y)  pu  3T  +  pv  3T  =  f(y)  u_  djo 
Tx  3y  Cp  3x 

+  1  J 3 |B(y)r^  KZ  3_  (CpT)l  +  jj_  f  3_u  1P 

fe(y)  rjCp  9y  1  cP  ay  J  cP|_3y  J 

+  1  r^Cp  3  1  B(y)r^ (-Cppv’T1 )  I  -  pu * v 1  3u 

BTy)  3y  I  .  J  Cp  3y 

+  1  pu3  CD  D^y')  ( B.  9 ) 

2B(y)  Cp  ^2 


( B3 )  Transformation  to  Levy- Lees  Variables 

The  transformation  of  Probstein-Elliot-t  and  Levy-Lees 
(Ref  31)  is  used  to  transform  equations  (B.5)>  (B.7),  and 
(B.9).  The  Levy-Lees  variables  £ , q  are  defined  as  follows: 

(x)  =  /x  pe  ue  pe  ro^dx  (B.10) 


and 


n(x,y)  =  pe  ue  ro^  t**  £_  dy 

2£  °  pe 


(3.11  ) 


The  partial-derivative  operator  in  the  new  coordinate  system 
(£,p)  becomes 


K 


ue  ue  rc 
21 


9— 

pe 


(B.12) 


(B.13) 
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Cvl Vfrfofr 


'wy ty wj rty%  r;  v\  r, v.  *r4  *r. r.  r.  *V«"iY*vrvvT*£ 


Equations  (B.10)  and  (B.11)  gives  a  transformation 

2  1 

=  pe  ue  ye  ro  0 
3x 


9 n  =  t^  ue  ro^i 


The  dependent  variable  F,  0,  and  V  are  defined  as: 


F  =  u  ,  9  =  T 

ue  Te 


pe  ue  ue  ro 


f(y)  Fnx  +  rJtJov 


Continuity  Equation 

Equation  (B.5)  may  be  written  as: 

3  {r^pv}  -  - f ( y )  3_  (r^pu) 

3y  3x 

Integrating  both  sides  from  wall  (y=0)  to  source  point  y 
in  the  field: 

rJpv  =  -  Iy  f ( y )  3__  (rJpu)dy  +  roJ  pto  vu 
0  3x 


Using  Equation  ( B . 1 3 ) 


l  ~ 

rdpv  = 


(5*(k)n  +  "*  I?)  </  f(y)(r°'l0U)ynan  +  r0 


Expanding, 


r°pv  =  -C 


k  (°/n  f 


n  f(y)(ro^pu)  2g  dp 

t^  ue  ro^p 

/n  f (y) (roJpu)  2j  dp\ 

t^uero^p  J 


+  roJ  pw  vio 


Using  Equation  (B.16) 

rJpv  =  -E  S_  /  /n  f  (y )  dn 


x  35  1° 


4yn  f(y)  rJE'in'j 
|„  3_/Vn  f(y) 

3u  V 


<2n  \  +  re15  pco  vu) 


This  equation  can  be  re-written  as: 
...i .  f*  /\  ^  r ^  r 


r"  pv  = 


f(y)  Fdr;) 

+  ro^  pqj  vu)  2 £  t^  -  2£  n  /n  f(v)  F'1' 

£  £  X  C 
sx 


(B.17) 


Simplifying,  the  first  two  terms  on  the  RHS  cf  “he  above 
equation  are  defined  as: 

v .  r. Q/n  f(y)^Fdn  +  roJ’  pco  -  Fv'  1 


Substituting  (B.18)  in  Equation  (3.17) 

i-jpv  -  £x  . rr  „  ,, 


JZT  t1 


v  -fe.  nx  Q/n  f(y)  Ed' 

^x 


(B.18) 


•^pvtJ*  =  v  -  2£  n  n/n  f(y)  Far 


v  is  given  by: 


^ r 

rJptJ  I  2£ 


2£  nx  f(y)  t 


(B.19) 
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solving  for  the  dependent  variable  V, 

V  =  2C  {rjpvtJ‘  +  2C  n  f(y)  F} 
E 


substituting  for  ?x  and  re-arranging, 

V  =  2g  (“f (y )  Fnv  +  rjtjov 


pe  ue  pe  ro 


2J 


“ 


2C 


which  is  the  same  as  Equation  (B.16). 


© 


To  obtain  the  continuity  equation,  differentiate  equation 

(B.17)  w.r.t.  n,  as  follows: 

Vn  =  -  2?  3_  {  2C  f (y)  F}  +  0 
35 


Simplifying 


Vn  =  -2?  f (y  -  f(y)  F 


Rearranging  the  above  equation  as 

Vn  +  f(y)F  +  2?  f (y)  =  0 


( B. 20 ) 


i  sdb 


The  continuity  equation  in  Levy-Lees  variable  including  the 
blockage  effect  is  obtained. 

Momentum  Equation 

Rewriting  the  momentum  equation  (B.7)  as: 


f(y)pu  3u  +  pv  3ui  =  -f(y)d£  +  j_  1  3 

3x  3y  dx  j  B (y )  dy 


(• 


i9u  -  pu'v' } 


B(y)r^ 


1  pu2  cd  D(y} 

2  B(y)  £2 


(B.21  ) 


Equation  (B.21)  is  divided  into  small  terms  I,  II,  III,  and 

IV  for  convenience.  Each  term  will  be  tackled  individually. 
Considering  term  I, 


I  -  f(y)  pu  3u.  +  pv  3u 


3x 

3y 

=  f(y)  pue 

F 

*x  It 

+  nx 

+  _P 

35 

3n 

ro^  pt ^ 

v  -  2£ 

nx 

f  (y)F 

nx  f (y) 

3  ue  F 

2C 

3n 

=  f(y)  pue2 

FCx  F{  + 

f(y)  pue' 

2  Fnx  Fn 

+  f (y)  Pue  F  ue^  +  p  t £ uero^  p 

roJptJ  2£ 

hlL  ue  F  -  2?  n  f(y)  ue  F  F 

^  r|  x  r) 

2£ 

I  -  f (y)  pue  F  £  F  +  f(y)  pue2/Fn  F 

s  x  n 

+  f(y)  pueF2Cx  ue^  +  ue2p  £xl  F^  -  f(y)nx  pue/FF 

2c  2e  n 

I  =  f(y)  pue  F  £x  F^  +  f (y )  pue  £x  F2ue^ 

2 

,  pue  E  v  _  , 

+  —  x  F^  (3.22) 

2£ 

Considering  term  II, 

II  =  -f(y)  d£ 
dx 


From  invicid  flow 


dpe  =  -pe  ue  due 
d£  d£ 
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Substituting  in  II  above, 


II  =  +f(y)  pe  ue  ue  =  f(y)  pe  ue  £  ue 

A  A  c, 


II  =  f(y)  pe  ue  £x  ue^ 


(B.23) 


Considering  term  III  now. 


III  =  1  1_ 

BTyT  rj 


_^B(y)  ,3u  -  pu'vj^ 


But ; 


eT  =  -p  u  ■'  >■  * 
3u/  3y 


Where  T  is  the  intermittency  factor. 


=  1 

1 

n  1 B(y)  r" 

(u  +  er) 

3u 

BTyT 

J 

y  I 

9y 

=  1 

1 

Hv  B(y)  r j 

( u  +  e  r ) 

ue? 

bTJT 

r*5 

y  L 

-  i 

1 

It^uero^p  1 

B  (  v )  r  j 

(u 

BTyT 

rj 

L  *  j 

n  n 


uet°  uerc 
2C 


'-1  pF^ 

i 


III  =  1 _  ue'to^o  {t^  (p  +  er)  o  F  } 

BTyT  or  "  ' 


n  "V  r%  r:  ^  v*r;  r;p;^r;iV^i 


I  '"■  ” 'W  '>  »T  ■ 


Finally,  considering  term  IV, 


IV  =  “I  cn  D^v)  P 
2  B(yU2 


■1  cn  D(y)  p  ue2  F2 
2  B(yH2 


2  t-,2 


IV  =  -1_  CD  D ( y )  p  ue'1  F 


B(yH' 


(B. 25) 


Substitute  Equation  (B.22)  to  (B.25)  in  Equation  (B.21) 

2 

f(y)  pue2  F  ?x  F?  +  f(y)  pue  ?x  F2  ue?  +  pve  gxV  Fn 

2? 


© 


f(y)  pe  ue  £x  ue?  +  ue3ro2jp  { B (y )  i2^  (y  +  er)oF„} 


B(y)  2C 


n  n 


-1  D(v )  pue2  F2 

2  B(yH2 


(B. 26) 


Simplifying : 


ue_  pF  £  {2£  f(y)F.  +  f(y)  2£  F  ueJ 

2£  x  ^  ^ 


ue 


A.  r 

+  pUe  lx  VF^  =  f(y)  pe  ue  ue.  +  ue3ro2j‘p 

2C  n  X  ?  BT7T2  f 


{t2j‘  (u  +  er)  pF^ } 


1  °n  p(y )  pue2F2 

2  B(yH2 


3  /•:> 
?1  sv 
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5 


i*  •**  */  O 


Dividing  the  above  equation  by  £  ue 

X  2£  p 

2K  f (y)  F  Fr  +  f (y)  2£  F2  ue.  +  VF 
*  ue  s.  h 


=  f(y)  ££  2^  uer  +  uero2^  fB(y)  t2^  (u  +  er)pF  } 
p  ue  *  B(y)?x  n 


-  1  "n  D(v)  M  F2 
2  B(y)t2  Cx 


Substituting  8,  as  defined  earlier  by  eqn  (35): 

2£  f(y)  F  F?  +  f(y)  BF2  +  VFn  =  f(y)6B 

+  uero2^p _  {t2^  B(y)(u  +  eT)pF  } 

2  i  MM 

B(y)peueuerc  J 


-  CD  D(.v)  2( 

B(y)t  peueuero^ 


Simplifying  further. 


2C  f (y)  F  F^  +  f(y)B(F2  -  9)  +  VFn  s 


__j _ _  |B(y)t2^(u  +  eT)  £_  1  F  1 

B(y )  pe  ue  ji 


-  cD  DixL 


al 


B(y)£  peueuero' 


=  1,  { B (y ) t2j*  ( 1  +  erU  Fn)  -  CD  D(y) 


bTTT 


zr 


B  ( y ) S,  peueuero 


2j 


'3  ■:*  t*  v 


T*  ^T^V'1 


1  {B(y)t2je  a  F  }  -  CD  D(y)  _ 


Fiyy 


B(y)£^  peueuero11 


Hence,  the  momentuin  equations  become: 


2K  f(y)FF-  +  F(y)B(F2  -  0)  +  VFn  =  1 _  { B(y) t2-Ui F  > 

t  n  ^yy  n  n 


2  9  * 

B(y)£,  peueyero~^ 


The  source-sink  term  is  defined  as  follows: 


4>ss  =  CD  l 

B(y)  £2  ro^peueye 


(3.23) 


To  obtain  the  final  form  of  mean-momentum  equation  including 


the  effect  of  surface  roughness,  in  the  Levy-Lees  variables, 


substitute  (B.28)  in  the  momentum  eqn  above, 


2K  f(y)FF.  +  f(y)B(F2  -  0)  +  VFn  =  1  {B(y)t2^e  A?  > 

K  n  bTTT  n  r 


-3>ssF' 


(B.29) 


Energy  Equation: 


Rewriting  the  Energy  equation  (B9)  as: 


f(y)pu  3T_  +  pv  3T_  =  f(y)  u^_  djc  +  1  .  9 

9x  3y  Cp  dx  B  ( y )  r^1  Cp3y 


[B(y)rJ’  KA  3_(CpT)  +  li/3u\2  +  1  ._3 _ 

CP  3y  J  Cp  ^ 3y  J  B(y )rJCp3y 


{ B ( y)r^  ( -Cppv '  T  ’  ) }  -  pu' v'  +  1  p u C D  D ( y ) 

Cp  3y  2B(y)  l2  Cc 


(3.30) 


1 


4 


II  =  f  (y)  u  d£  +  1  8  {B(y)rJ  U  3_  (CpT)}+  ^  /  Du  \ 

3x  ^THCpay  Cp  3y  Cp  y  3y  J 


But: 


=  K  Pe  u  u  a^d  Ka  =  jj_ 
dx  x  e  e  ec  Cp  pr 


Hence,  term  II  becomes: 


II  •■=  -f(y)  e  U  pe  »e  u  )  +  _J -  n  3 

CP  5  B(y)rJCp  ^ 


l'B(y)rJ'  jj_  Op  Te  n  6  }  +(u_ 


p?  e  y  n  'fe  eri  i" 


(B.32) 


III  =1  4  9  {B(y)r*^(-Cp  pv'T  r  )}  -  pu  ’  v '  3u 

B(y) rJCp  3y  Cp  3y 


But:  -Cp  pv'T1  =  +Kt  3T 


•pu'v'  =  c  c-U 


nv  3_  fB(y)rj(KT  3?))  +  c_  (du' 

B(y)rJCp  3n  3>’  CP 


■X_  |_  {B(y)  r3  K  Te  n  6  >  4  e  (ueF  n  ): 
3n  J  Cp 


B(y)rJCp 


(B. 33) 


IV  =  1  pu3  CD  D(y)  =  1  CD  D(^  pv  F3 


1  pu 

2BT7T 


2BT7T  0  2  Cp' 


( B .  3  4 ) 


‘  »  '  .  •  k'«  ,*'*>  . *•  N  .V  A  A 


Substitute  (B.31)  to  ( B. 34 )  in  IB. 30) 


f(y)pue  F  Te  £x  05  ♦  f(y)pue  £x  F6  T  *  pue  T^ve,, 

5  2£ 


=  -f(y)  u  F  (£  p  u  u  )  +  ~  . 

sr  x  e  e  ^  i^^i(B(y)rVpVy%> 


+  k-  (».  Fn  n„)2  +  2 


e  n  y 


%  3_(B(y)rJ  K  T  n  8  ) 

BTyTr^Cp  3n  I  e  y  n 


+  e  (u  F  n  )2  +  1  CD  3  -n  3 

C?  e  n  y  myT  mo2  °u€ 


Multiplying  the  above  equn  by  2 £; /pT  u  £  , 

G  v  A 


2K  f(y)  F  9,  +  25  f(y)  F_  6T  +  V< 

S  m  “  ** 


-f(y)  2CF  (p  u  u  )  +  nv  d_  B (y )r^  u_  Cp  ny(2c) 

T  Cp  e  e  e£  n/  x  i_  3n  Pr  puQ  C  n 

p  e  ^  -1 _ _ _ _ 

^  A 


+  2£u_  Ue  Fn  ny2  +  ny 
•  CP  P  ^  t>  ( ,,  v.  j  n. 


.?_[  3(j 
3n 


j  kt  2£  n 


v)rJ  T 


pu„  S, 


From  invicid  flow 


T 

o 


T  + 


e 


1 

1 


u 


e 


K 


or 


2£ _  T  =  _  2£  u 

CpTe  e?  CpTe  ue  e? 


_  0  =  -«e  (B.36) 

CpTe 


2?  f(y)F  +  V0n  =  IA  +  IIA  +  IB  +  IIB  +  IIIA  (B.37) 


XA  +  TB 


B(y) 


i _ a_  r  ? 

3yL 


(y)rJ  pi _ 


nv2(2£)  +  kt  2?nv2 


Pr  u  £, 
e  sx 


Cp  pue£x 


I  y 

i"  en 

3  X  J 


1  3__r  B(y)r^  pi_  + 

B(y)r*5  8n  L  Pp  °P 


Km  2Cn 


Xr-  e 


oue5x  n 


] 


But  =  e 
Cp  Prt 


1  3 

B  ( y )  r  ^  /  jj__  +  e 

\  _ 

2£  ro^  t2^' 

pu6 

2  1 

'  A  1 

B(y)rJ*  8n 

L  \?T  Prty 

'  2£ 

2  i 

pu  ro  J  p 
e  e 

'*  e 

%  nJ 

=  1 

3  T  t2^  B(y)rj  / 

^1  + 

e\?r  pu 

1 

e  1 

B(y) 

r0  3n  [  \ 

7Prt  'Ve 

Pr 

M 

=  1 

4 

t2irj  B(y)  e 

i 

el 

B(y)rJ* 

3n[ 

Pr 

nJ 
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Where 


£  =  1  +  e  Pr 


U  Prt 


£  =  ou 


P  u 

e 


—  2  2  2  2 
1IA  +  Up  =  2C  ^  ^e_  !n_V  +  2i£  u  V  n.y 


CP  pTe 


C'P  P  Te 


=  2£ 


ue  ny  ("  jj _  + 

PTe  ^x  LCP  CPJ 


2  i  .  2  i  p  2  2  _ 

2£  e  _ _ e  1  +  c_  I 

L  p  J 


PCPT  oC/"“2j 


e  2^(ro  J  pe  ue  ue) 


,  2j  U  ^  r— , _  2 

=  t  J  PU  0  { £ )  :• 


p  u  Or'’71 
w  e  e  •  e 


=  t^S,  a  e  F 


where 


e  =  1  +  e_ 
U 


and 


a  =  o 


CpT 


2  „3 


III*  =  1  °D  2£  PUe  F 


2BT7T  £2 


Re-arranging 


III, 


Cp  D(y)  uR' 


B(y H5  cPTe  ro2J 


EF 


p„  u  u 
e  e  e 
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Appendix  C 


Difference  Relations 


Three  point  central  differencing  relations  for  n  and 
three  point  or  two  point  backward  differencing  relations 
are  used  to  reduce  the  transformed  continuity,  momentum 
and  energy  equations  to  finite  difference  form.  It  is 
assumed  that  all  data  are  known  at  the  solution  station 
i-2  and  i-1  (Fig  3),  and  the  unknown  quantities  need  to 
be  determined  at  the  grid  points  for  "i"  station.  The 
finite  difference  equations  were  derived  with  the  stipu¬ 
lation  that  a  function  could  be  described  at  any  point  by 
i.  Taylor  series  expansion  about  another  point.  For  the 
present  study  the  approximation  was  made  that  for  any 
functional  value  F, 


F  ( i  —  1  ,  j  )  =  F  ( i ,  j  )  -  A£.  1  3F(i,j)  +  i-1  32F(i,,j)  (C.  la) 


F(i  +  1,j)  =  F(i,  j  )  -  (A£5_2  +  AC._i  )  3F(i,.1 ) 

3^ 


+  (ACi_i  +  ACi_2)"  9^F(i,.i)  +  - 

2 


(C. lb) 
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The  scaling  factor  for  £  derivative  are  defined  as  follows: 


Y1  ---  2(A£._2  +  2a£._1) 


(A£i_1  +  A?i_2) 


Y2  =  2(A£j__1  +  A£i_2) 


Tae-^T 


Y3  =  2A£' 


4?i-2(«i-i1+  A5i-2^ 


Y4  =  (Aei-2  *  AEt_i ) 
«i-2 


Y5  -  A£  . 


i-1 


A£ 


-2 


(G.2) 


(C.3) 


(C.4) 


(C.5) 


(C.6) 


Using  scaling  factors  defined  in  eqn  (C.2)  to  (C.6),  Equation 
( C . 1  a )  and  (C.lb)  can  be  solved  to  yield 


rail  =  Y3  F(i-2,.1)  -  Y2  F(i-1  ,j)  j  II  F(i,.j)  (0 

L  35  J  i. j  2A?i-1 


and 

F ( i , j )  =  Y4  F  ( i  - 1  ,  j  )  -  Y5  F(i-2,j)  (C.8) 


1  08 


'  © 


Terms  of  the  order  .  4?.  0  , ,  ,  , 

i-I  1-2  or  smaller  are  neglected. 

This  procedure  produces  truncation  errors  of  the  order  of 
^i-2  ^i-1  * 

2 

3F  3  F 

For  obtaining  expression  for  -7r~f  — w*  etc,  the  Taylo] 

dn  9p 

series  expansions  are  next  written  about  the  unknown  grid 
point  (i,j)  in  the  n-direction  as  follows: 


F ( i , j  + 1  )  =  F(i,i)  +  An,  F  (i,j)  +  An.i  F  ( i , j  )  +  ...  (0.9a) 


J  n 


A  2 

An  , 


F(  i  »  j  -1  )  =  F(i,i-1)  -  An.  1  Fn  ( i  ,  j  - 1  )  +  i-1  F  ( i  ,  j  - 1  ) .  . 

J “  '  p  MM 

(C . 9b) 


Equations  (ASa)  and  (A8b)  can  be  solved  to  yield  expres¬ 


sion  for  F  and  ?  as  follows: 

nn  n 


F  ( i »  j  )  =  Y6  F  (  i  ,  j  +1  )  -  2Y7  F(i,.i)  +  Y8  F(i,.1-1) 

1 1 1 1  *  P  p 

An  j  Ay  A^  j 


(C.10) 


yi.j) 


=  Y9  ?(i,.i+1)  -  Y10  F  ( i  ,  .j  )  -  Y  8  F(i,.j-1  ) 


2  An  . 
0 


2An  • 

J 


(C.11  ) 


where  the  coefficients  Y6,  Y7.....Y10  coefficients  are  defined 
as  follows: 


r-  AP|  -i  1 

2/fl  + 

J 


( C  .  1  2  ) 


vWuVa'  av‘a'  '.'  V  i'  OV\  •  -  \ 


(C.  1 3 


Y7  =  An  ,/An  .  1 

J  J 

Y8  =  2/  An.1-1ri  + 

An,,  L 


An , 


Ar>i 


(C.  1 4 


Y9  =  2/ 


[1 


_^L_- 

Anj-1- 


(C.15 


Y10 


1 


An, 


An 


( G .  1 6 


2 

Finite -difference  representations  for  _36_,  9  6  ,  98 ,  3V  are 

9n  «  2  3£  9n 

derived  simultaneously.  dn 

All  of  these  are  substituted  in  equations  (30)  to  (32) 
to  obtain  the  corresponding  finite-difference  equation.  Du 
to  their  recurring  appearance,  the  following  symbols  and 
their  definitions  are  introduced: 


FM1  =  Y4  F( i-1 , j )  -  Y5  F( i-2 , j ) 

TM1  =  Y4  T ( i-t , j  )  -  Y5  T(i-2, j  ) 

and  ( C . 1 7 

FM2  =  Y2  F ( i-1 , j  )  -  Y 3  F(i-2,j) 

TM2  =  Y2  T (i-1  ,  j )  -  Y3  T(i-2, j  ) 


1 1  0 


Non-linear  terms  of  the  form, 


), 


where  G  and  H  represent 


any  typical  dependent  variable,  appearing  in  the  governing 


equations  need  to  be  linearized  in  order  to  obtain  a  system 
of  linear  difference  equations.  Quantities  of  this  type 
are  linearized  by  using  equation  (C.7)  and  (C.8),  i.e., 


(G  ffh.J  =  (YA  Gd-M)  -  Y5(i-2,j)HY3  H(i-2,j) 

-Y2  H(i-I.J)  +  Y4  H(l,j)}/2A5i_1  +  C  AC._2)  (C.18) 

For  example. 


(F  ||)i>j  =  FK1  { Y2  TV'  -  T5  FM2}  (C.1 9) 

The  procedure  used  to  linearize  these  r. :n-linearized  product 

terms  such  as  /_9G\ /  is  as  follows: 

^  9n  M  3n  ’ 


-  (iii) 

1  v  3n  n  3n  •  ,  ^  3n  • 


^  U 
( ' 


-  (  9H  \ 

n  •  '  3r.  •  ,  .  ,  V  3n  J  .  ,  . 


1,J  "■  i-1  ,j  d’'  n  .  n  '  ^  ’  -'-I  <  9n  , 


'i  /  3  G  \ 

v  3n  1  ,  1  3n  ; .  . 

1  - 1  * i » j 


(  C  .  20  ) 


where  the  terms  / 3G\ 

’h-i.i 


and  ^  3H  s 

Sn  i-1 ,! 


are  discretized 


according  to  eqn  (C.11)  but  evaluated  a~  proceeding  station 


111 


(i-1).  The  linearization  for  quantities  of  the  form 
is  obtained  by  setting  H  to  G  in  eqn  (C.19). 


2 


For  example: 


3F(i 


=  2FY  i 


-  FY 


,1.21) 


F2  =  2?(i , j )  F ( i  —  1 , j )  -  F ( i— 1 » j )  (C.22) 


where  FY  denotes 


and  F(i,j)  are  unknown. 


a  knovm  quantity  whereas 


All  terms  had  been  represented  in  finite  difference 


form,  and  the  final  step  incorporated  these  linearized  models 
into  equations  (32),  (31 ) ,  and  (32)  to  derive  the  overall 
system  of  finite  difference  equations  (Ref  19:67-71). 


Appendix  D 


Coefficient  in  Difference  Equation 

Equation  (53)  to  (55)  are  the  difference  equations  used 
to  represent  the  partial  differential  equation  for  the  con¬ 
servation  of  momentum  and  energy  and  continuity  respectively 
These  equations  are  repeated  here  fo?’  convenience. 

Momentum : 


A1  (n ,  1  )  Fn-1  +  A1(n,2)  Fn  +  A1(n,3)  Fn  +  1 

+  B1(n,1)  Gn-1  +  B1(n,2)  6n  +  Bl(n,3)  0n+1 
+  Cl  (n ,  1  )  Vn-1  +  C1(n,2)  Vn  +  C1(n,3)  Vn  +  1  =  D1(n)  (B.1 


Energy: 


A2 (n , 1 )  Fn-1  +  A2(n,2)  Fn  +  A2(n,D  Fn+1 

+  B2 (n , 1 )  0n-1  +  32(n,?'  0n  +  B2(n,3)  0n+1 
+  C2(n,1)  Vn-1  +  C2(n,2;  'n  +  C2(n,3)  Vn+1  =  D2(n)  \L,2 


Continuity : 


A3(n,-')  Fn-1  +  A3(n,2)  Fn  +  A3(n,3)  Fn  +  1 

+  B3(n,1)  0n-1  +  B3(n,2)  0n  +  B3(n,3)  0n+1 
+  C3(n,1)  Vn-1  +  C3(nf2)  Vn  +  C3(n,3)  Vn+1  -  D3(n)  (B.3 


These  equations  are  obtained  fror.  equations  (48),  (49),  and 
(50)  and  the  difference  quotients  presented  in  Appendix  'C'. 
The  coefficients  A1(n,1),  Bl(n,1)  and  so  forth  in  equations 
(D.1),  (D.2),  and  (D.3)  are  functions  of  known  quantities 
evaluated  at  station  i-1 ,  and  i-2. 


The  coefficients  are  as  follows: 


A1  (N  ,  1  ) =Y8#XL* (2 .  *XLM1  WEM1  'DY-OlLMI  *YYM1  *IM1  #XLPM1  rTY+ 

BYPP*EM1  *XLM1  /BYM1  -VM1  ) ) _ 

A1  (N  ,  2  )  =  -  ( 4.  ^XLttXLM1  ,>fEM1  ^Y7/DY+2  .  “XL"  (  XI'!'  “BY Ml  +BYPP" 
EMI  *XLM1  /BYM1  +EM1  #XLPM1  #TY-VM1  )  #Y1  0-2  .  ••BX2“FM1  -  ( XFY  ( N ) 

XBE+XSS  (N)  )  ^SEPtSEP^  ( 2 . *Y 1  "FM1  -F:'2  )  "jlFY  (Y  N  -X  ) _ 

A 1  ( :N ,  3 )  =  X  L *(  2  .  *  X  L  M 1  *  E  ff 1 ~ Y  6  7  D  Y  +  ( X  I,' '  1  “  I Y  Ml  - 1 ? '  1  *  X  LP  M 1  *  T  Y  + 

BYPP*EM1«XLM1/BYM1-VM1 )«Y9) _ 

~B1  (N,1  )=-XL»EM1*XLPM1*FYtfY8 _ 

~5i~{N,2;=DX2*XBE»XFY(N)-2.  *XLPM1  *FY;"Y1  0  _ 

B1  (N ,  3  )  =XL~:fEM1  *XLPM1  *FY*Y9 _ 

Cl (N.1 )=C1  (N,3)  =  0. 


Cl  (NJ[_2)=DX2*FY 


[  =  -2.  *XL‘;tXAL'"XLM1  “EMI  -FY-Y8 


A2(N,2)  =  -(4.*XL*XAL*XLM1  *EM1  •;-'FYi  Y1C  +  B 
XFY  ( N  )-3.**nX2ttXAL»XSS(N)*?M1 


TO- 


.  Y1  WTM1  -TM2  ) ' 


■  -  «ct: 


A2  (H ,  37=2 .  '"XL^XAL':fXLM1  *EM1  "?Y: 

B2  ( N ,  1  ) =XL#Y8  “  ( 2  .  #XLM1  -ETM1  /  ( F3;' BY )  -  ~ 
XLPM1 *ETM1  «TY+BYPP»XLM1  “ETMI  /  BY Ml  - ? 3 
B2  (n  ,  2 )  -  ( 4 .  "XL"XLM1  *ETM1  /  (F3. 

XLPM1  #ETM1  *TY  +  BYPP‘:tXLM1  #ETM1  /  BY 


•byT+ 

mi  -?: 


'BY  Ml  +2.- 

?H) _ 


•7 


::ymi  +2.* 

•XL“Y1  0*2.0/ 


PR+SEP*X»Y1*XFY(N)*FM1 ) _ 

B2(N,3)=XL*(2.  #XLM1  *ETM1  *Y6/DY+  ,XLM1  rrZ?YMT  -2.  *XL?M1  “ 

:  ETM1  *TY+BYPP*XLM1  *ETM1  /BYM1  -PR-'Y"1  )  "Y9  '?“ 

C2(N,2)=-TX2*TY  _ 

C2(N,1 )=C2(N.3)=0. 

~A  3  ( N ,  1  )=A3(-T,  3J  =  0. _ 

A3  N,2)=(DX2+X*Y1  )*XFY(N) _ 

B3  (N  ,  1  )  -B3  ( N ,  2  f=B3l N ,  3  )  =  0  .  _ 

C3(N,1  )=-XL*Y8  _ 

C3(N,2)=-2.*XL*Y1  0 
C3(N. j)=XL»Y9 

di  TTFmwmm  *xlpmi  *ty-vmi  )  - f i  “ “ . xb i-'-xfy  U)*ox2 

+X*YHDX2*XSS  (N )  )*SEP  _ 

D2(n)  =  DX2 *(XLPM1  *ETM1  #TY/?P.-VM1  ;  T Y  +  : X 2  ' X A L '”;X LM 1  *EM1  * 
TY##2-X*XFY(N)  *Y1  *TM1  *FM1  #XFY(K )  *SE?  +  :>X2*XAL*2  .  XSS(N)  * 


FM1  *^3*SEP _ 

D3(N)=X*FM2*XFY(N)' 


(P.4) 

(D.5) 

(D.6) 

(D.7) 

(D.8) 

(D.9) 

(D.10) 

(D.11  ) 

(D.12) 

(P.13) 

(D.14) 

(P.1n) 


(P.16) 

(P.17) 

7  D .  1  «  1 
(D.19) 
(^.20) 
(D.21  ) 
(P.22) 
(P.23) 
(P.24) 
(P.25) 

(D.26) 

(P.27) 

(D.28) 
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k.  _  a  —  » 


The  additional  quantities  appearing  in  equations  ( D . 4 )  to 
(D.28)  are  defined  as  follows: 


5 


FM1  = 

Y4 

F(i-  1 

-  Y5 

F(i-2, 

■  3 ) 

(0.29) 

TM1  = 

Y4 

T(i-1 

*.i ) 

-  Y5 

T  (i-2 , 

■  3 ) 

(D.30) 

VM1  = 

Y4 

T  ( i-1 

.3) 

-  Y5 

V ( i-2 , 

»3 ) 

( D.  31  ) 

FM2  = 

Y  2 

F  ( i-1 

..1 ) 

-  Y3 

F( i-2 , 

-  j  ) 

(D.32) 

TM2  = 

Y2 

T(i-1 

» j ) 

-  Y3 

T  ( i-2  , 

•  j  ) 

(D.33) 

(pu  )/  (pu  ) 

1 

e  1 

tfhich 

may 

be  w] 

fitten 

as 

XLM1  =  TM1 


1  +  ( 


S> 


'e  l 


TM1  +  VT 

e  l 


(air  only)  (D.34) 


I  <f> 


t 

S’ 

A 


XLPM1  =  £p  =  XLM1 
2TM1 


(f-) 

e  i 


(f-> 

1  e  i 


+  TM1 


(1.35) 


EMI  =  e(i,j)  =  s  ( i  ,  ,i  - 1  )  +  e  ( 1 , ;;  ~  +  f.  ( i ,  ,j  + 1  )  (D.36) 


rr 

i 

\  « 

i 


W 


ETM1  =  e  ( i  >  j  ) 


BYM1  =  B(y) (i,j )  (eqn  20) 


FY  =  _Y9  F(i-1,j+1)  -  Y1_0  F  (  i  - 1  ,  j  ) 
2An  An 


Y 8  F(i-1 , j-1 ) 
2  An 


(D.37) 

(0.38) 

(0.39) 


i 
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-“r-3  -  -  •  •  '  •  .  V  -V-" 

— V  _  M  _  k  k  - - M  -  m  -  »  x. 


TY 


EYM1 


3TYM1 


BYPP 


Il_  T  ( i  —  1  ,  j  +  1 )  -  Y1_0  F(  i-1  ,  j  )  -  Y8_  F(i-1,j 
2  An  An  2  An 


Y9  e(i-1,j+1)  -  TMO  F(i-1,j)  -  Y8  e  (i-1  .  j 
2  An  An  2An 


I2_  e  (i-1  .  j  +1  )  -  Y10  e  ( i-1  ,  j  )  -  Y8_  c  (i-1  #  j 
2An  An  2 An 


Y9  BY  (i-1  ,.1-H  )-Y10  BY  ( i-1  ,.j  )-Y8  EY(i-1,j-1 

2An  An  2An 


XEE  =  8  ( i , j )  =  ,2£  due \ 

^ u  d£  ' 


XAL  =  o ( i , j ) 


e 


i 


XI.  =  Lt 


XFY  =  f (y ) ( i , j )  ( eqn  21) 


xss 


s  s  ( i ,  j  ) 


(eqn  37) 


1)  (0.40) 

1)  (D. 41 ) 

1)  (D.42) 

)  (D. 43 ) 

(0.44) 

(0.45) 

(0.46) 

(0.47) 

(D. 48) 


Equations  (0.1),  (0.2),  and  (D.  3 )  are  solved  simultaneously 
by  Gaussian  elimination  technique. 


AeB#ndi><_E 


PROGRAM  LISTING 


PROGRAM  I TRACT (INPUT, OUTPUT , TAPES- INPUT , TAPE6-0UTPUT) 

c  ********************************************************** 

C  THIS  PROGRAM  SOLVES  COMPRESSIBLE  TURBULENT  BOUNDARY-LAYER 

C  INCLUDING  THE  EFFECT  OF  SURFACE  ROUGHNESS, THE  SHAPE. HEIGHT 

C  AND  DENCITY  OF  THE  ROUGHNESS  ELEMENTS  IS  SPECIFIED  TO  THE 
C  CODE  IN  SUBROUTINE  RUFVAR. 

C  ********************************************************** 

COMMON  G,  PR,  REY,  XMINF,  OMEGA,  BO,  TW,  P10,  TlO,  RIO.  VISIO,  TE. 

1  PE.  RE, UE, VISINF, SU, EPS, DS, DYW,  SI , ERROR, TC, TA, I EDGE, IEND1 , INTACT , 

2  PRT, XXK, BTRX, XLAM, VARPRT , X INTER, SEPO, XCHS(B) , IPRN (9) . EQ (300) , 

3  EN(300) , EP (300) , ETO (300) , ETN(300) , ETP (300) , FO (300) , FN (300) , J2DA, 

4  FP(300) , TN(300) , TO (300) , XNN (300) , VN(300) , VO (300) , VP (300) , TP (300) . 

5  Di (300) ,D2(300) ,D3(300) 

DIMENSION  Y (300) , A1 (300, 3) , A2 (300. 3) , A3 (300, 3) , B1 (300.3) . 

1  B2 < 300 , 3 ) , B3 ( 300 , 3 )  Cl (300, 3) , C2 (300, 3) . C3 (300, 3) 

COMMON/CPDATA/  CP (24) , XP (24) , DP (24) , IPRES 
COMMON/BLRVAR/BYO (300) , BYN(300) , BYP (300) . XFY(300) , XSS(300) . 

1  XBRE, XHRE, XLS, XCD 

1100  FORMAT < 1H0, 17X, 3HS/L, 15X.2HCP,  15X,  6HP/PINF) 

1101  FORMAT ( IX, 3 (4Y, E15. 9) ) 

200B  FORMAT (IX, *PROFILE  FAILED  TO  RELAX  AT  M  -  *,I5) 

8002  FORMAT  (3E10.6) 

S003  FORMAT  (1015) 

9002  FORMAT (1H1.47X* INTERACTING  BOUNDARY  LAYER  SOLUTION*) 

9003  FORMAT  (7H0GAMMA-F6. 3, 4H  PR-F6.3.5H  MFS-F6.3.7H  REYFS-E10. 4 . BH  TFS 
1 (K) -F7. 1 , 1 1H  B0-TW/T10-F6.4.5H  EPS-FB.5) 

9004  FORMAT (5H0P10-, E10. 4, 7H  RH010«,E10. 4,5H  T10-, E10. 4 . 7H  VIS10-.E10.4 
1, 4H  SI-, E10. 4) 

9005  FORMAT (7H00MEGA-, F7. 4, 2X, 6HPRT  «  , F7. 4, 2X , 7HBTRX  -  .F7.4) 

9019  FORMAT (10X,*WITH  INTERMITTENCY  CORRECTION*) 

9020  FORMAT < 10X, *WITHOUT  INTERMITTENCY  CORRECTION*) 

9021  FORMAT (10X, ITWO-DIMENSIONAL  BOUNDARY  LAYER*) 

9022  FORMAT ( 10X, *AXISYMETRICAL  BOUNDARY  LAYER*) 

9023  FORMAT <  10X.  /* - SOLUTION  NEEDS  MORE  POINTS  FOR  CONVERGENC— */) 

C 

C  INPUT  INITIAL  CONDITIONS 
C 

READ (5, 0002)  G,  PR,  XMINF,  REY,  TA 
READ (5, 8002)  DS,  SI,  OMEGA,  ERROR.  XXK 
READ (5, 8002)  BO, BTRX,  PRT,  X INTER,  DYW 

READ <3, 8003)  IEDGE,  INTACT,  IDIFF,  IEND1 ,  MSP,  J2DA  .  IPRES 
READ (3, BOOS)  (ICHS(I),  I  -  1,  B> 


ono  non  non 


READ (5, 8003)  (IPRN(I),  I  -  1 ,  9) 

XLAM«.3*BTRX 

IF(IPRES.EQ.O)  GO  TO  20 

READ (3, 8002)  DPMAX 

READ (5, 8002)  (CP(IJ) , IJ-1, IPRES) 

READ (5, 8002)  (XP(IJ) , IJ-1, IPRES) 

WRITE (6, 1100) 

XMSQ-XMINFIXMINF 
DO  10  IJ-1, IPRES 
PDPINF-1 . 0+0. 5*G*XMSQ*CP <  I J ) 

WRITE (6, 1101)  XP(IJ) ,CP(IJ)  .PDPINF 
10  CP ( I J) -PDPINF 

CALL  BMTHPR ( BTRX , DPMAX , G , XMSQ > 

COMPUTE  NOND I MENS I ONAL I Z I NG  QUANTITIES 

20  Zl-  1.  +  (G  -  1. )/2.*XMINFS*2 

P10  «  (1./<G*XMINF**2))*<Z1**(G/(G-1.))) 

T10  »  ( 1 . / ( <G  -  1. )*XMINF**2))*Z1 
RIO  -  G*P10/(T10*(G  -  1.) ) 

TINF  «  T10/Z1 
TW  «  B0*T10 

IF (OMEGA  .EQ.  0.)  GO  TO  101 
ViSlO  “  T 10* IOMEGA 

EPS  *  <<<G  -  1.  )*XMINF**2)** (OMEGA/2. )) /SQRT(REV) 

*'1SINF  -  TINF**OMEGA 
GO  TO  102 

101  'C-198. 6/  ( (G-l . ) *XMINF**2*TA) 

VIS10  -  <T10m.5)*<l.  +  TC)  /  (T10+TC) 

FPS  *  ( ( ( (1 .  +  ( 198. 6/TA) ) * ( ( (G  -  1 . ) *XMINF**2) **1 . 5) ) / < (SG 1 . >  *X 
JMINF**2)+(198.6/TA) ) )/REY)**.5 
VISINF  -  (TINF**i.5)*(l.  +  TC) / (TINF+TC) 

102  SU-198.& 

OUTPUT  INITIAL  CONDITIONS 
WRITE (6, 9002) 

WRITE (6, 9003)  G,  PR,  XMINF,  REY.  TA,  BO,  EPS 
WRITE <6, 9004)  P10.  RIO,  T10,  VIS10,  SI 
WRITE (6, 9005)  OMEGA, PRT, BTRX 
IF (XINTER. EQ. 1 . )  WRITE (6,9019) 

IF ( X INTER. EQ, 0. )  WRITE (6 , 9020) 

IF (J2DA. EQ. 0)  WRITE (6, 9021 ) 

IF < J2DA. NE. 0)  WRITE (fe, 9022) 

INPUT  INITIAL  PROFILE 

12  MSTART-2 

C  INITIALIZE  THE  STREAMWISE  LOCATION 
S-SI 

DS2-DS1-DS 

DX2DS-DX1DS-DXDS-0. 


1 1  7 


n  n  cj  non  non  nnn 


SEPO-1. 

C  INITIALIZE  THE  STREAMWISE  LOCATION 
Y  < 1 ) "0. 0 
DO  201  LL-2,300 
DY-XXK**(LL-2)*DYW 
201  Y (LL) -Y (LL-l ) +DY 
DO  700  LL-l,  300 
D1 (LL)-D2(LL)-D3<LL)-XNN<LL)-0. 

VP (LL) -VN  (LL)  -VO  (L.L)  —  Y  (LL) 

FP(LL) -FO(LL) "FN(LL) -TP (LL) -TN(LL) -TO (LL) -EP (LL) -EQ(LL) "EN (LL) - 
1  E7P (LL) -ETO (LL) -ETN (LL) -BYP (LL) -BVO (LL) -BYN (LL) -1 . 0 

700  CONTINUE 

DO  701  J  -  1,  300 
DO  701  I  -  1,  3 

701  A1 (J, I)»A2(J, I)-A3(J,I)-B1 (J,  I)-B2(J, I)-B3(J, I) -Cl (J, I) 

1  -C2(J. I)-C3(J, I)»0. 

PREF«G*XMINF**2 
TREF  «  (G  -  l.)*XMINF**2 

INITIALIZE  COUNTERS 

ICOUN-MSTART 
IQ-IEDGE 
1 0-1 
t  r  =  i 
lN^CH-0 
ITCNT1  «  1 
I IN-0 


*■**■  BEGIN  FIRST-ORDER  TRIDIAGONAL  MATRIX  SOLUTION 


f>  *  * 


DO  115  M-MSTART, IEND1 

IF (M. EQ. MSTART)  MP-MSTART 

IF (M.EQ. IEND1 )  MP-M 

IF(M. EG. (M/MSP) *MSP>  MP-M 

S-S+DS2 

DX2DS  »  DX1DS 

DX IDS  =  DXDS 


COMPUTE  LOCAL  PRESSURE  AND  PRESSURE  GRADIENT 

CALL  PRESSM  <S, XMINF , G, PBG1 , DPBG1 , TETNF, XME) 

COMPUTE  LOCAL  EDGE  PROPERTIES 

PE  -  PBG1/PREF 

P*  -  DPBG1/PREF 

TE  -  TETNF/TREF 

UE  -  SORT <2. 1 (T10  -  TE)) 

RE«G*PE/((G-1.0)*TE) 

TR-SU/ (TETNF*TA) 

IF  (OMEGA)  642,676,642 
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n  n  o  n  o  no  o  non  n  n  o  n  n  oon 


642  XNUE-TEMOMEGA 
GQT0688 

676  XNUE-TEM1 . 3*  <  1 .  +198. 6/ (TAtTREF) )  /  (TE+198. 6/ <TA*TREF) ) 

688  CONTINUE 

COMPUTE  LOCAL  XI  AND  STEP  LENGTHS 

DXDS-RE*UE*XNUE 
IF ( J2DA. NE, 0)  DXDS-DXDS*S**2 
IF  <M. EQ. 2)  DX1DB»DX2DS-DXDS 

DX2-. 34DS2* ( ( 1 . +DS2/D81 ) 4DX1DS+DS1 *DXDS/ (DS1+DS2) -DS2IDS2IDX2DS/ 

1  (DS1*  <DS1+DS2) ) ) 

REYNDE«RE*UE*S/XNUE 
REYEXT-REYIVI3INF4REYNDE 
IF (M.EQ. 2)  DX1-DX2 
IF (M. EQ. 2)  X«DXDS*SI 
X-X+DX2 

COMPUTE  STEP  LENGTH  FUNCTIONS 

Yl*2.  XDX1+2.  *DX2)  /  (DX1+DX2) 

IF < IDIFF  .EQ.  1)  Y 1  «  2. 

Y2* ( (DX1+DX2) /DX1 ) 12. 0 
Y3“(DX2*DX2/ (DXlt (DX1+DX2) ))*2.0 
Y4*>  (DX1+DX2)  /DX1 
Y3-DX2/DX1 
TWTE  -  TW/TE 

COMPUTE  ALPHA,  BETA,  AND  LAMBDA 

DUEDX— PP/(RE*UE*DXDS> 

XAL«UE*UE/TE 
XBE-2.0*X*DUEDX/UE 

ASSIGN  ROUGHNESS  VARIABLES  TO  ELMATX 
CALL  RUFVAR(X, XNUE) 

6998  LENGTH=IEDGE 

ASSIGN  THE  MATRIX  ELEMENTS  FOR  THE  FINITE  DIFFERENCE  EQUATIONS 
CALL  ELMATX (  M, DX2, X, XAL, XBE, TR, IDIFF, Y1 , Y2, Y3, Y4, Y5, TWTE, ITCNT1 . 
1  A1,A2,A3,B1,B2,83,C1,C2,C3) 

ASSIGN  THE  MATRIX  ELEMENTS  FOR  THE  FINITE  DIFFERENCE  EQUATIONS 


MATRIX  INVERSION,  SOLVE  FOR  F,  THETEA  AND  V 

CALL  MATEQN3 (FP, TP, VP, D1 , D2, D3, A1 , B1 , Cl , A2, B2, C2, A3, B3, C3, 3, LENGTH 
1  ,300) 

C 

C  MATRIX  INVERSION,  SOLVE  FOR  F,  THETEA  AND  V 


noon  on 


ITCNT1-ITCNT1+1 

N-IEDGE+1 

DY-DYWIXXKM  (IEDGE-2) 

VK-(VP(IEDGE)/(XXK*(1.-H.  /XXK) > -VP < IEDGE-1 >*< 1 . -1 . /XXK) IXXK- 
1  VP(IE0GE-2)*XXK/(1.+1./XXK))/DY 
DY-DYW*XXK**< IEDGE-1) 

KON-N+3 
DQ63IQ«N,KQN 
DY-DYW*XXK**(N-1)+DY 
FP(IQ)-TP<ia>-1.0 
65  VP (IQ) -VP < IEDGE-1 >+vK*DY 

C  INITIATION  OF  SIMILAR  SOLUTIONS 

IF<M.EQ. 2)  GO  TO  8020 
GO  TO  801 8 

8020  DO  8019  I-t,K0N 
VO ( I ) -VN ( I ) -VP ( I ) 

FO ( I ) -FN ( I > -FP  < I ) 

8019  TO ( I ) -TN ( I ) -TP ( I ) 

C  INITIATION  OF  SIMILAR  SOLUTIONS 

8018  IQ-IEPGE+1 

U  AND  THETA  PROFILES  ITERATIONS 
TAU2- (FP(2) -FP <  1 ) ) /DYW 
IF (ITCNT1 ,EQ. 2)  TAU1-10. *TAU2 
RT12-TAU1/TAU2-1. 

TAU1-TAU2 

IF  < ITCNT 1  .LE.  100)  GO  TO  7005 
WRITE (6. 2008)  M 
CALL  EXIT 

7005  IF (ABS (RT12) =GT; ERROR)  GO  TO  699S 
U  AND  THETA  PROFILES  ITERATIONS 


COMPUTE  BLT,  BDT (DELTA  STAR)  AND  BMT (THETA/ 

35  CO-TP (1) 

TPI*0. 

BLT-BLDT«BLMT^O. 

XNN ( 1 ) -0# 

DO  57  N-2, KON 
CY-DYW*XXK*t (N-2) 

C-TP(N) 

TPI“TPI+. 5KDY* (CO+C) 

CO-C 

XNN(N>-TPI*SQRT<2.*X)/<RE*UE> 

IF  < J2DA. NE. 0)  XNN (N> -XNN  <N) /S 

BLDT -BLDT + ( 2 . -FP ( N )  /TP(N) -FP(N-1 ) /TP(N-l) ) K XNN (N) -XNN (N-l ) ) /2. 
BLMT-BLMT+ (FP  <N> * ( 1 . -FP (N> ) /TP  <N) +FP  <N-1 ) * ( 1 . -FP (N-l ) > /TP (N- 1 ) ) 

1  *(XNN(N)-XNN(N-l))/2. 

IF (BLT.GT.O. )  GO  TO  57 

IF (FP (N) . GE.0.99)  BLT-XNN (N) - <FP(N> 99) * (XNN (N) -XNN (N-l > ) 

1  / (FP (N) -FP (N-l ) ) 

57  CONTINUE 


."V 


BLT«BLT*EP8 

BLDT-BLDTIEPS 

BLMT-BLMTtEPS 

C  COMPUTE  BLT,  BDT( DELTA  STAR)  AND  BMT (THETA) 

C 

C  COMPUTE  THE  EDDY  VISCOSITY  CUEFFICIENT 
IF(S.LE.BTRX)  GO  TO  56 

CALL  REYSTP  (KON, TR, X, TREF, XNUE, XBE, S, ITCNT1 ) 

C  COMPUTE  THE  EDDY  VISCOSITY  COEFFICIENT 

C 

58  ITCNT1*1 
C 

C  ASSESMENT  OF  GRID  PON ITS  IN  ETA 

C 

IF (INDCH)  71,  71,  732 

71  CONTINUE 

IF <M  -  20)  732,  732,  72 

72  IF(ABS (FP ( IEDGE-15) -FP ( IEDGE- 16) ) -0.0001)  73,73,74 

73  IF (ABS (TP ( IEDGE-15)  -  TP < IEDGE-16) >  -  .0001)  732,  732,  74 

c  t***tmmim*mmtm*tm*mmmmm*m******* 

C  FOR  TURBULENT  FLOW  SOLUTION  NEEDS  TO  BE  MONITORED  CLOSE 
C  TO  THE  OUTER  EDGE. WHILKE  ADDING  MORE  POINTS  THE  MAXIMUM 

C  VALUE  OF  I EDGE  BE  RESTRICTED  TO  MAX.  DIMENSIONS  MINUS  25. 

c  *********************************************************** 

74  IF (IEDGE. GT . 275)  GO  TO  75 
IEDGE=IEDGE+1 

IQ-IQ+1 
WRITE <6. 9023) 

75  DY»DYW*XXK** ( I EDGE-2) 

Y< IEDGE)  -  Y ( IEDGE-1 )  +  DY 
732  IQ  -  IQ  -  1 

C  ASSESMENT  OF  GRID  PONITS  IN  ETA 

C 

C 

C  COMPUTE  WALL  STRESS  AND  HEAT  TRANSFER  AND  OUTPUT  STATION 

CALL  CFSTNO  (TR, XNUE, X, S, XBE, M, BLDT, BLMT, BLT, PBG1 , DPBG1 , REYEXT. 
1  XME, MP) 

C  COMPUTE  WALL  STRESS  AND  HEAT  TRANSFER  AND  OUTPUT  STATION 

C 

C 

C  SHIFT  PROFILES  BACK  ONE  XI  STATION 

C 

NN  «  IQ  +  5 
DO  118  N*1,NN 
FN(N)«FO(N) 

FO(N)-FP(N) 

TN(N)«TO(N) 

TO  (N)  *TP  (N) 

VN (N) =V0 (N) 

VO(N)«VP(N) 

ETNvM)»ETO(N> 

ETO(N) «ETP(N) 


h 


•1 


v'i 


.1 


fel 

■V 


■8 
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( ICOUN,  IP,  IG,  IQ,  MSTART,  I  IN,  M,  S,  Y,  BLT,  XME) 


■■  .V 

v.-' 


EN<N>-E0(N) 

EO(N)*EP(N) 

BYN(N)»BY0(N) 

BYO<N)-BYP<N) 

118  CONTINUE 
DX1-DX2 
081*082 

IF  <M+1-ICHS( IG) )  114,113,114 

113  082*2. 0IDS1 
IG  -  IG+1 
INDCH  -  1 

IF  (M.EO.IENDl)  GO  TO  237 
GO  TO  111 

114  082*081 
INDCH  *  0 

IF  (M.EO.IENDl)  GO  TO  237 
GO  TO  111 
237  1 IN  *  1 
111  CALL  PRNCHS 
115  CONTINUE 
STOP 
END 

SUBROUTINE  PRESSM <8, XM, G, P, DPDX, T, YM) 

COMMON/CPDATA/  CP (24) , XP(24> ,DP(24> , IPRES 
100  FORMAT (  5X,*WARNING.... CALCULATION  IS  OUTSIDE  OF  THE  PRESCRIBED  PR 
1ESSURE  DATA.  S  IS  LESS  THAN  XP(i>*> 

200  FORMAT (  5X. CWARNING. ...  CALCULATION  IS  OUTSIDE  OF  THE  PRESCRIBED  PR 
1ESSURE  DATA.  S  IS  GREATER  THAN  XP(END)*) 

300  FORMAT (1X,3E15.9) 

IR-0 

IPM1-IPRES-1 
IF ( IPRES. EQ.O)  GO  TO  40 
DO  20  1*1, IPRES 
IF <8.LT.XP(1) )  WRITE<6. 100) 

IF(S.GT. XP( IPRES) )  WRITE (6, 200) 

IF (8.LE. XP<1) )  IR-1 
IF(IR.NE.O)  GO  TO  30 
IF<8.GE.XP(IPM1) )  IR-IPRES 
IF  < IR.NE. 0)  GO  TO  30 

IF < (8. GE.  XP ( I ) ) . AND. <S. LT. XP <1+1 ) )  >  IR-I 
IF (IR. EQ.O)  GO  TO  20 
:  SEEKING  THE  BEST  FIT 

RS* (S-XP(I ) ) / (XP( 1+1 ) -XP( I ) ) 

IF (RS.GT.0.5)  IR-I+1 
:  SEEKING  THE  BEST  FIT 

IF(IR.NE.O)  GO  TO  30 
20  CONTINUE 

30  IF (IR.GT. IPM1 )  IR-IPM1 

IRP-IR+1 
IRM-IR-1 

IF(IR.EQ.l)  IRM-IR+2 

:  COMPUTE  THE  CUBIC  SPLINE  COEFFICENTS 
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X1-(XP(IRP)+XP(IRM)-2.0*XP<IR))*<XP<IRP)-XP<IRM)) 

X2= (XP (IRM) ~XP  (IR))*(XP(IRM)  -XP  <  IR)  ) 

X3*(XP (IRP) -XP( IR) )  t  (XP ( IRP) -XP ( IR) ) 

X4*XP  < IRM) -XP  < IRP) 

X5*XP ( IRM) -XP ( IR) 

X6-XP  < IRP) -XP(IR) 

DETS«X5*X6*X4 

C2- (DP ( IR) tXl+DP (IRP) HX2-DP ( IRM) *X3> /DETS 
C3* (DP  < IR) <X4-DP < IRP) KX5+PP  < IRM) 1X6) /DETS 
C  COMPUTE  THE  CUBIC  SPLINE  COEFFI CENTS 
DXP-S-XP(IR) 

DXP2«DXP**2 

DXPF-DXP/20. 

DPDXl-DP(IR) 

P*CP(IR) 

DO  10  1*1,20 

X-ItDXPF 

X2«X*X 

DPDX2-DP ( IR) +C2*X+C3*X2 
P«P+0. 5* (DPDX1+DPDX2) *DXPF 
10  DPDX1-DPDX2 

DPDX*DP  < IR) +C2*DXP+C3*DXP2 
T«P**( (6-1.0) /G) 

YM«SQRT<2.0*( <2. 0+<G-l .  0)  >XM*XM)  /  (2.  0*T> -1 . 0)  /  (G-l .  0) ) 

WRITE (6,300)  S,P« DPDX, T, YM 
GO  TO  50 
40  P-1.0 

DPDX-O. 

T=P**< (G-1.0)/G> 

YM-SQRT (2. 0* ( (2. 0+ (G-l . 0) *XM*XM) / (2. OtT) -1 .  0)  /  (G-l .  0) ) 

DO  RETURN 

END 

SUBROUTINE  SMTHPR (BTRX, DPMAX, G, XMSQ) 

COMMQN/CPDATA/  CP (24) , XP (24) , DP (24) , IPRES 
100  FORMAT (IX, *FIRST  CP  DATA  POINT  YIELDS  ADVERSE  PRESSURE  GRADIENT  TO 
10  STEEP  FOR  CALCULATION  TO  CONTINUE*) 

200  FORMAT ( 1H0, 1 1 X, 3HS/L, 15X, 2HCP, 1 1 X, 6HP/PINF, 14X, 4HDPDX) 

300  FORMAT  <1X,4(4X,E15.9) ) 

DPT0L-DPMAXU.01 

C  COMPUTE  THE  TRAILING  EDGE  DPDX 
IPM1«IPRES-1 
IPM2-IPRES-2 
DX1-XP ( IPM1 ) -XP (IPRES) 

DX2-XP  < IPM2) -XP< IPRES) 

DX12-DXUDX1 

DX22-DX2IDX2 

DP ( IPRES) - (CP < IPM2) IDX12-CP ( IPM1 ) *DX22-CP (IPRES) *  <DX 12-DX22) ) / 

1  (DX1IDX2* (DX1-DX2) ) 

C  COMPUTE  THE  TRAIL  NG  EDGE  DPDX 
10  IMAX-0 

C  COMPUTE  THE  LEADING  EDGE  DPDX 
DX1-XP  <2) -XP ( 1 > 
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DX2=XP(3) -XP  < 1 ) 

DX12-DX1IDX1 

DX22-DX24DX2 

DP (1)* (CP (3) IDX12-CP (2) *DX22-CP< 1 >  % (DX12-DX22) ) / (DX1 >DX2> (DX1- 
1  DX2) ) 

IF (DP ( 1 ) . GT> DPMAX )  WRITE (6, 100) 

IF (DP(1) . GT.DPMAX)  CALL  EXIT 
C  COMPUTE  THE  LEADING  EDGE  DPDX 
DO  20  I*2,IPM1 
IM1-I-1 
IP1-I+1 

DXl-XP(IMl)-XPd) 

DX2-XP(IPl)-XPd> 

DX12=DX1*DX1 

DX22-DX24DX2 

DP  <I)»  <CP(IP1 ) *DX12-CP( IM1 ) IDX22-CP ( I ) * (DX12-DX22) ) / (DX1 *DX2* 

1  (DX1-DX2) ) 

20  IF ( <DP(I)  .GT.DPTOL)  .AND.  (XPd) . LE.BTRX) )  IMAX-I 
IF(IliAX.EQ.O)  GO  TO  50 

C  SMOOTHING  THE  CP  DATA  IN  THE  LEADING  EDGE  REGION 
IMM1*IMAX-1 
IMP1-IMAX+1 
DX1°XP(IMM1)-XP(IMAX> 

DX2*XP ( IMP1 ) -XP (IMAX ) 

DX12-DXUDX1 

DX22-DX2CDX2 

CP ( IMM1 ) * (CP ( IMP1 ) *DX12-CP( IMAX ) * (DX12-DX22) -DX1 *DX2* (DX1-DX2) 

1  < DPMAX) /DX22 
GO  TO  10 

C  SMOOTHING  THE  CP  DATA  IN  THE  LEADING  EDGE  REGION 
50  WRITE (6, 200) 

DO  30  1*1 , IPRES 

PC*2. 0* (CP (I) -1 . 0) / (G*XMSQ) 

30  WRITE (6, 300)  XP (I > , PC, CP (I ) , DP< I ) 

RETURN 

END 

C 

c  ***************  *************** 

c 

C  THIS  S/R  HAS  BEEN  MODIFIED  TO  INCLUDE  THE 

C  EFFECT  OF  SURFACE  ROUGHNESS  999999999999 

C 

c  *3*************  **************  **************** 

SUBROUTINE  CFSTNO  (TR, XNUE,  X ,  S,  XBE, M, BLDT, BLMT, BLT. PBG1 , DPBG1 . 

1  REYEXT,XME,MP) 

COMMON  G,  PR,  REY,  XMINF,  OMEGA,  BO,  TW,  P10,  TiO,  RIO,  VIS10,  TE, 

1  PE, RE, UE, VISINF, SU, EPS, Dw. DYW,  SI , ERROR, TC, TA, 1EDGE, IEND1 , INTACT, 

2  PRT , XXK, BTRX, XLAM, VARPRT , XINTER, SEPO, !CHS(8) , IPRN (9) , EO (300) , 

3  EN (300) , EP (300) , ETO (300) , ETN (300) ,  ETP (300) , FO (300) , FN (300) , J2DA, 

4  FP(300) ,TN(300) , T0(300) , XNN(300) , VN(300) , VO (300) , VP (300) , TP (300) , 

5  Dl (300) ,D2(300) ,D3(300) 

COMMON/BLRVAR/BYO (300) , BYN (300) ,  BYP(300) , XFY(300) , XSS(300) , 
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1  XBRE, XHRE, XLS, XCD 

2000  FORMAT < 1H0, 10X, 5HS/L  -,E15,8> 

2001  FORMAT  <21 , 7HXME  «,  E15. 8* 2X,  7HPE  E15. 8, 2X, 7HDPPINF«,E15.8, 

1  2X? 7HX8F  ", E1S.8, 2X, 7HTW/TE  -,E15.B) 

2002  FORMAT <2 <f7HBLT  E15. 8, 2X, 7HBLMT  », E1S. 8, 2X, 7HBLDT  «,E15.8, 

1  2X, 7KPEYMT  ■, E15. 8, 2X , 7HREYDT  -,E15.B> 

2003  FORMAT <2X,7HCFN0  «, E15. 8, 2X, 7HCFEN0  «, E15. 8, 2X, 7HSTN0  -,E15.B, 

1  2X, 7HSTEN0  ■, E15. 8, 2X, 7HREYEXT-, E15. 8) 

TAURbO. 

FPSI-O. 

FP(1>«0. 

C0-FP(l)t*2 
DO  100  N»2,IEDGE 
XYY«XNN(N)*EPS 
IF (XYY-XHRE)853,853.854 

853  DY»DYW*XXK**(N-2> 

C«FP(N)**2 

FPSI«FP8I+.5*DY*(C0+C) 

TAUR-FPSI*.3*XCD*XBREtUE*SQRT (2. *X> /XLS**2 
CO«C 

100  CONTINUE 

854  TWTE*TP(1> 

IF (OMEGA. EQ. 0. )  GO  TO  855 
IF (OMEGA  .EQ.  1.)  GO  TO  8551 
XLM1  =  1. / (TWTEtf ( 1 .  -  OMEGA)) 

GO  TO  856 

855  XLM1  =  ( .  O+TR )  tSQRT  ( TWTE )  /  ( TWTE+TR ) 

GO  TO  856 

8551  XLM1  «  1. 

856  CONTINUE 

Yll-(<2.+XXK)*(l.+XXK+XXK**2)+l.+XXK>/< ( 1 . +XXK) * < 1 . +XXK+XXK**2) ) 

Y12«(1.+XXK+XXK**2)/XXK**2 

Y13- (1 . +XXK+XXK**2> / <XXK**3* ( 1 . +XXK) > 

Y14«1./(XXK**3*(1.+XXK+XXK**2>  > 

TAU-XLMl*R£*XNUE*UE*UE*(-Yil*FP<l>+Y12*FP<2)-Y13*FP(3)+Y14fFP<4)  ) 
1  / (DYWtSQRT (2. IX) ) +TAUR 

QS  »  XLM1*RE*XNUE*UE*TE*<Y11*TP(1)-Y12*TP<2)+Y13*TP(3)-Y14*TP(4) ) 
1  / (DYWtSQRT (2. *X) IPR) 

IF ( J2DA. NE. 0)  TAU»TAU*S 
IF < J2DA. NE. 0)  QS=QS*S 
STNO  *  0. 

IF (BO  .NE.  1.)  STNO  =  EPSIQS/((1.  -  BO)*(TE  +  .5*UE**2)) 

STENO  ■  STNO/ <RE*UE) 

CFNO  -  2. tEPSITAU 
CFENO  «  CFNO/ (REt JE*UE) 

REYDT-REYEXTIBLDT/S 
REYMT-REVEXT*BLMT/S 
C  SELECTION  OF  THE  OUTPUT 
IF(M.NE.MP)  GO  TO  1000 
SELECTION  OF  THE  OUTPUT 

OUTPUT  STATION  DATA 
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WRITE <6, 2000)  S 

WR I TE  < 6 , 200 1 )  XME , PBG 1 , DPBG 1 , XBE , TWTE 
WRITE <6, 2002)  BLT, BLMT , BLDT, REYHT, REYDT 
WRITE <6, 2003)  CFNO, CFENQ, STNO, STENO, RE YE XT 
1000  RETURN 
END 

SUBROUTINE  ELMATX  (  M,DX2, X, XAL, XBE, TR, IDIFF, Y1 f Y2,  Y3, Y4,  Y5,  TWTE, 

1  ITCNT1,  A1 „ A2,A3,B1,B2,B3,C1,C2,C3> 

COMMON  G,  PR,  REY,  XMINF,  OMEGA,  80,  TW,  P10,  T10,  RIO,  VIS10,  TE, 

1  PE, RE, HE,  VISINF,  SU,  EPS,  DS,  DYW,  SI ,  ERROR,  TC,  TA,  I  EDGE,  IEND1,  INTACT, 

2  PRT,  XXK,BTRX,  XLAM, VARPRT, X INTER, SEPO, ICHS(B) , IPRN(9) , EG (300) , 

3  EN (300) , EP (300) , ETO (3C0) , ETN (300) , ETP (300) , FO (300) ,  FN (300) , J2DA, 

4  FP ( 300 ) ,  TN ( 300 ) , TO ( 300  > , XNN ( 300 ) , VN ( 300 ) , VO ( 300 ) ,  VP ( 300 ) ,  TP ( 300 )  , 

5  D1 (300) , D2 (300) , D3 (300) 

CQMMON/BLRVAR/BYO (300) , BYN (300) , BYP (300) , XFY (300) ,  XSS (300) , 

1  XBRE, XHRE, XLS, XCD 

DIMENSION  A1 (300,  3) , A2 (300, 3) , A3 (300, 3) , Bi (300, 3) , B2 (300, 3) , 

1  B3  (300, 3)  ,  Cl  (300,3)  ,  C2  (30C,  3) ,  C3  (300, 3) 


:  THE  INNER  EDGE  BOUNDARY  CONDITION 

% 

DO  8011  1-1,3 

801 1  A1  ( 1 , 1 )  *A2  (1,1)  «A3  <  1 ,  I ) *B1  ( i  ,  I)  *B2 (1,1)  «B3 (1,1) *>C  1(1,1)  =*C2  (1,1) 
1  »C3 ( 1 , I ) *0. 

Ai < 1 , 1 ) *1 . 0 
B2<1.1)*1.0 
D1 < 1 ) *0. 

D2 ( 1 ) *TWTE 

IF  (SEPO.  EG).  0.  )  GO  TO  8012 

C3 ( 1 , 1 ) *1 • 0 

D3(l)«0. 

GO  TO  8013 

8012  XL-DX2/(2.0*DYW) 

A3(l, 1>-DX2+X*Y1 

C3( 1 , 1 ) *-2. *XLt(2.+XXK)/(l.+XXK> 

C3(1,2)«2.»XL*(1.+XXK)/XXK 
C3 ( 1 , 3) *-2.  *XL/(XXK» (l.+XXK) ) 

D3 ( 1 ) *0. 

* 

:  THE  INNER  EDGE  BOUNDARY  CONDITION 

:  THE  FIELD  POINTS  EVALUATION 


8013  NM1-IEDGE-1 

DO  8014  N«2, NM1 
DY-XXKt* (N-l ) *DYW 
DYM1-DY/XXK 
XL-DX2/ (2. 0*DY) 

Y6-2. / ( 1 . +DYM1 /DY) 

Y7-DY/DYM1 

Y8"2. / ( (DYM1/DY) ♦ ( 1 , +DYM1 /DY) ) 
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Y9=2. / ( 1 . +DY/DYM1 ) 

Y10-1.-DY/DYM1 
SEP- 1.0 

IF(FO(N).LE.O. )  SEP-O. 

IF  < ITCNT1  .GT.  1)  GO  TO  7000 
IFUDIFF  . EQ.  1)  GO  TO  7501 
FMI-Y4IF0 (N) -Y5tFN (N) 

TM 1 "Y4  $  TO ( N ) - Y5  >  TN ( N ) 

VMl-Y4tVQ(N) -YStVN (N) 

IF (SEPO. EQ. 0. )  VMl-VO(N) 

EMI  -<Y4lt<EQ<N-l>+E0(N>+E0<N+l))-Y5*<EN(N-l>+EN<N>+EN<N+l) )  )/3. 
ETM1  - < Y4t (ETQ (N-l ) +ETQ (N) +ETO (N+l >  > -Y5*  <ETN (N-l ) +ETN (N>  +ETN (N+l 
1  ))>/3. 

BYM1  - ( Y4* ( BYQ  <N- 1 >  +B YQ ( N  >  +B YD ( N+ 1 } > -Y5* ( BYN ( N- 1 ) +BYN « N ) +BYN ( N+ 1 
1  )))/3. 

GO  TO  7001 
7501  FM1  “  FO(N> 

TM1  -  TO(N) 

VM1  *  VO <N) 

EMI* (EO (N-l ) +E0 (N) +E0(N+1 ) ) /3. 

ETM1- (ET0(N-1 ) +ETQ (N) +ETQ(N+1 > > /3. 

BYM1- (BYO (N-l ) +BYO  <N) +BY0(N+1 ) ) /3. 

GO  TO  7001 

7000  FM1  *  FP (N) 

TM1  *  TP (N) 

VM1  *  VP (N) 

EMI* (EP (N-l ) +EP  <N) +EP (N+l ) > /3. 

ETMl«(ETP(N-l)+ETP(N>+ETP<N+l>>/3. 

BYM1- (BYP (N-l ) +BYP (N) +BYP (N+l ) ) /3. 

7001  IF (OMEGA  .EQ.  0.)  GO  TO  684 
IF (OMEGA  .EQ.  1. )  GO  TO  6841 
Xt,Ml*l./(TMltt(l. -OMEGA) ) 

XLPM1- (OMEGA-1. )*XLM1/TM1 
G0TO625 

6841  XLM1-1. 

XLPM1-0. 

60T0623 

684  XLM1- ( ( 1. +TR) ISQRT (TM1 ) / (TMl+TR) ) 

XLPMl-XLMlt (TR-TM1 ) / (2. *TM1 * (TMl+TR) > 

625  IF ( ITCNT1 . GT. 1 )  GO  TO  626 
FY*(Y9*FO(N+l)/2.-Y10tFO(N)-Y8*FO(N-l) /2. ) /DY 
TY- (Y9*T0(N+1) /2. -Y10*TO(N) -Y8KT0 (N-l ) /2. ) /DY 
EYM1* (Y94E0 (N+l ) /2. -Y1O4EM1-Y8IE0 (N-l ) /2. ) /DY 
ETYM1- ( Y9*ET0 (N+l ) /2. -Y10*ETM1-Y8*ET0 (N-l )/2.) /DY 
BYPP«(Y9*BY0<N+l)/2.-Y10fBYMl-YS*BY0<N-l>/2. >/DY 
GO  TO  627 

626  FY- ( Y?*FP (N+l ) /2. -Y10*FP <N> -Y8*FP (N-l ) /2. ) /DY 
TY*(YV*TP(N+1) /2.-Y10tTP(N) -YBITP(N-l) /2. ) /DY 
EYMl»(Y9*EP(N+l>/2.-Y10*EMl-YB*EP(N-l) /2. ) /DY 
ETYM1- (Y9IETP (N+l) /2. -YlOIETMl-YBtETP (N-l > /2. ) /DY 
BYPP- (Y9JBYP (N+l ) /2. -Y10*BYM1-YB*BYP (N-l ) /2. ) /DY 

627  IFdDIFF.EQ.  1)  GO  TO  7502 


o  o  o  n  o  o 


FM2=Y2*F0  <N) -V3IFN (N) 

TM2-Y2*T0(N>-Y3*TN<N> 

60  TO  7505 
7502  FM2  -2.*FQ<N) 

TM2  -2.*T0(N) 

7505  CONTINUE 

A1  <N,l>*Y8*XL*<2.*XLM*EMl/DY-<XLMl*EYm+EMl*XLPMl*TY+BYPP* 

1  EMlIXLMl/bYMl-VMl) > 

A1  <N,2>  —  <4.*XL*XLMi*EMl*Y7/DY+2.*XL*<XLMl*EYMl+BYPP*EM*XLMl 
1  /BYM1+EM1IXLPM1 tTY-VMl ) IY10+2. *DX2*FM1 * (XFY (N) IXBE+X93  <N) >  *SEP* 

1  SEP*  <2.  *Y1*FM1-FM2)  *XFY  <N) *X) 

A1 <N,  3) “XL* (2. *XLM1*EM1 IY6/DY+  <  XLM1 *EYM1+EM1 *XLPM1 tTY+BYPP* 

1  EMl*XLMi/BYMl-VMl)*Y9) 

B1 <N, 1>h-XL*EM< *XLPM1 IFYSYS 
B1 <N, 2)UDX2*XB2*XFY  <N)-2. *XL*EM1  *XLPKi *FY*Y10 
B1 <N,3)«XL*EM1*XLPMI*FY*Y9 
Cl (N, 1)-C1 (N,3)-0. 

Cl  (,M,2)  —  0X2IFY 

A2(N, l)*-2. *XL*XAL*XLM1*EM1*FY*YB 

A2(Nf2)«-<4.*XL*XAL*XLMl*EMt*FY*YlC+SEP*X*(Yl*TMl-TM2>* 

1  XFY(N)-3. *DX2*XAL*XSS  <N) *FM1**2. *SEP) 

A2  <N, 3) “2. *XL*XAL*XLM1*£M1*FY*Y9 

B2(N, 1 ) »XL*Y8* (2. *XLril*ETMl/ (PRIDY) - (XLM1 *ETYMl+2. *XLPM1*ETM1*TY 
1  +BYPP*XLMl*ETMl/BYMl-Pn*VMl) /PR) 

B2(N, 2) ■- (4. *XL*XLM1*ETM1*Y7/ (PR*DY) + ( XLM1 *ETYMl+2. *XLPM1 *ETM1 *TY 
1  +BYPP*XLM1*ETM1/BYM1-PR*VM1)*XL*Y10*2. 0/PR+SEP*X*Yl tXFY <N) *FM1 ) 
B2(N,3)»XL*(2.*XLMl*ET-ll*Y6/DY+(XL^11»ETYhi+2.*XLPMUETMl*TY 
1  +BYPPIXLM1 *ETM1 /BYMl-PRtVMl ) *Y9> /PR 
C2<N  2)— DX2*TY 
C7.  . 1)«C2(N?3>»0, 

A3(N, l)-A3<N,3)-0. 

A3(N,2)«(DX2+X*Y1)*XFY(N) 

B3(N,1>-B3(N,2)-B3<N,3)*0. 

C3<N,  D—  XL*Y8 

C3(N,2)«-2.*XL*Y10 

C3<N,3)*XL*Y9 

D1 (N) =DX2*FY* (EMI *XLPM1 *TY-VM1 ) -FM1 XX2*  ( XBEKXFY (N) *DX2+X*XFY (N) *Y1 
1  +XSS (N) *DX2) *SEP 

D2(N>*DX2* ^M1*ETM1*TY/PR-VM1)*TY+DX2*XAL*XLM1*EM1*FY**2-X*Y1 
Is...  -M'  ,  .  (N)*SEP+DX2*XAL*2.*XSS(N)  *FM1**3*SEP 
D3(N) ■X*rM2*XFY <N) 

8014  CONTINUE 

THE  FIELD  P^ts  EVALUATION 


THE  OUTER  ED6E  BOUNDARY  CONDITION 
DO  8015  1*1,3 

8015  A1 (I EDGE,  1)*A2 (IEDGE, I)  “A3  HEDGE, I)*B1 (IED6E, I)*B2(IEDGE, I)*B3( 
1  I EDGE , I ) -C 1 ( I EDGE , I ) *C2 ( I EDGE, I ) *C3 ( I  EDGE, I ) *0. 

A1  HEDGE,  3)  *1.0 


o  n  o 


B2(IEDGE,3>*1.0 

D1 (I EDGE) *1.0 

D2(IEL'GE)-1.0 

IF (SEPO. EQ. 0. )  GO  TO  8016 

XL-DX2/ <2. tDYWtXXK** (IEDGE-l ) ) 

FM2-Y24FO ( IEDGE) -Y3IFN ( I EDGE) 

IF(IDIFF.EQ.l)  FM2*2.*F0< IEDGE) 

A3(IEDGF,3)-DX2+X*Y1 

C3 (IEDGE, l)-2.IXXK*t3*XL/(l.+XXK> 

C3(IEDGE,2)  —  2,*XXK*(1.+XXK)*XL 

C3( IEDGE,  3) *2. IXXKIXLt (2. IXXK+1 . ) / (1. +XXK) 

D3(IEDGE>»X*FM2 
GO  TO  8017 

8016  VM1«V0( IEDGE) 

IF(ITGNT1.GT. 1)  VM1=VP (IEDGE) 

C3 ( IEDGE, 3) "1 . 0 

D3(IEDGE) *VM1 

8017  CONTINUE 
C 

C  THE  OUTER  EDGE  BOUNDARY  CONDITION 
RETURN 
END 

SUBROUT  I NE  PRNCHS  (  I COUN .  I P,  I G,  I Q ,  (1ST  ART ,  1 1 N ,  M ,  S ,  Y ,  BLT ,  XME ) 

COMMON  G,  PR,  REY,  XMINF,  OMEGA,  BO,  TW.  P10,  T10,  RIO,  VIS10,  TE, 

1  PE,RE,UE,VISINF,SU,EPS,DS,DYW,SI,ERROR,TC,TA, IEDGE, IEND1 , INTACT, 

2  PRT, XXK, BTRX, XLAM, VARPRT, X  INTER, SEPO, ICHS (8) , IPRN (9) , EO (300) , 

3  EN(300) , EP (300) , ET0(300) , ETNC300) , ETP (300) , FO (300) , FN (300) . J2DA, 

4  FP (300) ,  TN (300) , TO (300) , XNN (300) , ON (300) , 00(300) , VP (300) , TP (300) . 

5  D1 (300) , D2 (300) , D3 (300) 

DIMENSION  Y (300), Z (7, 16) 

25  FORMAT  ( 1 HO , 45 X , 23HPR0F I LE  FOR  STATION  S  -F14.8) 

40  FORMAT (8H0N*  15F8.4  ) 

41  FORMAT (8H  ETA*  13F8.4  ) 

42  FORMAT (BH  FI*  15F8.4  ) 

43  FORMAT (8H  Tl*  15F8.4  ) 

44  FORMAT (8H  VI*  15F8.2  > 

46  FORMAT (8H  EO*  13FB.2) 

507  FORMAT (8H  Y/BLT*  15F8.4  ) 

509  FORMAT <8H  R0/R0E*15F8. 4  ) 

510  FORMAT (8H  ML/ME*  15F8.4) 

311  FORMAT (8H  PT/PTE-15FB. 4  ) 

312  FORMAT <8H  PT/PE-  13FB.4  ) 

513  FORMAT <8H  H/HE-  15F8.4  ) 

IF(ICOUN-IPRN(IP) )  31,38,31 

OUTPUT  PROFILE  DATA 

38  KONT-IQ-1 
J2-0 

WRITE (6, 23)  S 
D050J1*1,KQNT, 15 
J2-J2+1 
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KQN=J2*15 

WRITE  (6,40)  (XNN(N) , N-Jl , KQN) 

WRITE  (6,41)  <Y(N),N-J1,K0N> 

WRITE  (6,42)  (F0(N) ,  N-J1,K0N> 

WRITE  (6,43)  <T0(N>,  N«J1,K0N> 

WRITE  (6,44)  ( VO ( N ) ,  N-J1,K0N> 

WRITE(6,46)  <EQ(N),N-J1,K0N) 

I-Jl-1 

IF (M.EQ.MSTART)  60  TO  50 
D0530JX-1, 15 
I-I  +  l 

Z<1, JX)-EPS*XNN(I)/BLT 
2(2, JX)-FO(I) 

C  2(3, JX)-TO(I) 

2(3, JX)-1./T0(I) 

PTPED-(G-1.0)*TE*T0(I) 

IF  (PTPED)  777,777,778 

777  PTPED- 1. 

778  Z  <4, JX)-UEf Z  <2, JX) / (PTPED) I*. 5 
PTPE-2 <4,JX)*Z(4,JX) 

IF <2 <4, JX) -1.0) 504, 504, 505 

504  PTPE* ( 1. 0+ ( ( (G-1.0) /2. 0) *PTPE) )$*(G/(G-1.0)> 

PTEPE-(1. 0+  ( ( (6-1 . 0)  /2.  )  *XME**2>  >  **  <6/  (6-1 .  >  > 

G0TQ506 

505  PTPE-  ( <  (6+1 . 0)  4PTPE/2. 0)  **  (6/  (G-l .  0) ) )  *  ( ( (G+i .  0)  /  ( (2. 0*G*PTF'E)  -  (G- 
11.0)))**<1.0/<G-1.0)>) 

PTEPE«<  ((6+1.  )*XME**2/2. )  *t  (6/  (6-1 . ) ) )  *  ( <  <6+1 . )  /  ( (2.  *G*XME**2S-<6- 
&1.  )))**<!. /(G-l.))) 

506  2  (5,  JX)  =PTF'E 

C  Z  <6, JX )  -  PTPE4PE/P10 
2 (6, JX) -PTPE/PTEPE 

2(7,  JX)*(TE*T0(I)/(UE*UE)+.5*F0(I>*F0(I)>/(TE/(UE*UE)+.5> 
2(4,JX)-Z(4,JX)/XHINF 
530  CONTINUE 

WRITE (6, 507)  <Z(1,N) ,N«1, 15) 

WRITE <6, 509)  <Z (3, N) , N-l , 15) 

WRITE (6,510)  <Z(4,N),N-1, 15) 

WR1‘TE<6, 51 1 )  (Z  <6,  N) ,  N-l ,  15) 

WRITE  <6, 512)  <  Z (5,  N) , N-l , 15) 

WRITE <6, 513)  (Z  <?,N> ,N*1, 15) 

50  CONTINUE 

IFdIN  .EQ.  1)  RETURN 
ICOUN-O 

51  ICOUN-ICOUN+1 

IF (M+l-ICHS (16) )  3601,3600,3601 

3600  IP- I P+1 
ICOUN-IPRN(IP) 

3601  CONTINUE 
RETURN 
END 

SUBROUTINE  REYSTR  (KON, TR, X, TREF, XNUE, XBE, S, ITCNTI ) 

COMMON  6,  PR,  REY,  XMINF,  OMEGA,  BO,  TW,  P10,  T10,  RIO,  VIS10,  TE, 


1  PE,  RE,  UE,  VISINF,  SU,  EPS,  DS,  DYW,  SI ,  ERROR,  TC, TA,  IEDGE,  IEND1 ,  INTACT. 

2  PRT,  XXK, BTRX,  XLAM, VARPRT, X INTER, SEPO, ICHS(B) , IPRN(9) ,ED(300) , 

3  EN  <  300 ) ,  EP  < 300 ) ,  ETO  <  300 ) , ETN ( 300 > , ETP ( 300 ) , FO  <  300 ) ,  FN ( 300 > ,  J 2DA , 

4  FP (300) ,  TN (300) ,  TO (300) , XNN (300) , VN(300) , 00(300) , VP (300) , TP  <300) , 
3  D1 (300), 02(300), 03(300) 

TTR« (TA+1 12. ) / (TAtTREF+l 12. ) 

CO-'P(l) 

DD-EP ( 1 ) -XNN ( 1 ) ■TPI-BLT-O. 

SHEAR  STRESS  AT  THE  WALL  AS  THE  SCALING  FUNCTION 

Yll-(<2.+XXK)*<l.+XXK+XXK**2)+l.+XXK>/< ( 1 . +XXK) t < 1 . +XXK+XXK**2) ) 

Y12«(1.+XXK+XXK**2)/XXK*<2 

Y13- ( 1 . +XXK+XXKIK2) / ( XXK**3t ( 1 . +  XXK) ) 

Y14«l./<XXK**3t<l.+XYK+XXK**2>> 

FETW* (-Y1 1*FP( 1 ) +Y12IFP (2) -Y13HFP (3) +Y14XFP (4) ) /DYW 
FETW=ABS(FETW) 

XLM1W* v ( 1 . +TR) ISQRT (TP (1 ) ) / (TP ( 1 ) +TR> > 

PI2-XLM1W*F£TW 

SHEAR  STRESS  AT  THE  WALL  AS  THE  SCALING  FUNCTION 

DO  1  N**2,  KON 

DY-DYW*XXK**(N-2) 

XLM1" ( < 1 .  +TR) ISORT (TP (N) ) / (TP (N>  +TR) ) 

C-TP(N) 

TPI«TPI+.3*DY*<C0+C) 

CO-C 

XNN (N) “TPI USQRT  (2.  *X)  /  (RE*UE) 

IF ( J2DA. NE. 0)  XNN <N) “XNN (N) /S 
IF (BLT. GT. 0. )  GO  TO  2 

IF (FP (N) . GE. 0. 99)  BLT=XNN(N) - (FP (N) 99) * (XNN (N) —XNN < N— 1 ) ) 

1  /(FP(N)-FP(N-l) ) 

DD-DD+<<i.-FF(N)>*TP<N>+(l.-FP<N-l) ) *TP(N-1) ) IDY/2. 
Pn«=S0RT(2.tX*REY/(TREF**1.5ITTR)  )*TPI**2/(XNUE*TP(N)  **3) 

IF ( J2DA. NE. 0)  PI1-  1/S 

0Y«DYW*XXK**(N-1 

DYM1-DY/XXK 

YS“2. / ( (DYM1 /DY) * ( 1 . +DYM1/DY) ) 

Y9*2. / ( 1 . +DY/DYM1 ) 

Y10*l . -DY/DYM1 

PI2-XL.M1IEP (N) *ABS  <Y9*FP (N+l ) /2. -Y10*FP (N) -Y8*FP (N-l ) /2. ) /DY 

CEB ICE-SMITH-MONSINK I 8  EDDY  VISCOSITY  MODEL 
YPLUS«SQRT (PI1IPI2) / (26. tXLMl) 

IF (YPLUS. GT. 50. )  YPLUS-50. 

EP  (N>«*.  16tPIl$  <1. -EXP (-YPLUS) >  t*2#ABS ( Y9* 

1  FP (N+l ) 12. -Y10JFP (N) -YS*FP (N-l ) /2. ) / (DY*XLM1 ) 

CEB I CE-8M I TH-MONS I NK I 8  EDDY  VISCOSITY  MODEL 

TRUNCATE  THE  INNER  REGION  CALCULATION 
IF (EP (N) . LE. EP (N-l ) )  EP (N) =EP (N-l ) 

TRUNCATE  THE  INNER  REGION  CALCULATION 

CONTINUE 
DO  3  N“1 . KON 


n  o  n  n 


XLMl-< (l.+TR)*SQRT(TP(N) ) /<TP(N)+TR> ) 

DD1«. 016BIS0RT (2.  *X*REY/ (TREF**1.3*TTR) )  *DD/ <XNUE*XLM1*TP(N) M2) 

IF<J2DA.NE.O)  DD1-DD1/S 

IF (DDl.LE.EP(N) )  EP(N)«DD1 

XXXX-. 412* ( (S-BTRX) /XLAM) **2 

IF(XXXX.GT.50. )XXXX-50. 

EP  <N)«EP(N> * ( 1 . -EXP (-XXXX) ) 

IF (XINTER.EQ.O. )  EP(N>*1.*EP(N) 

IFtXINTER.EQ.  1. >  EP(N)«i.+  ( 1 . 75/ ( 1 . +5. 5* (XNN (N) /BLT) **6) +1 . ) * 

1  EP (N)  /2. 73 

3  ETP  <N) *1 . +PR* (EP (N) -1 . ) /PRT 
C 

RETURN 

END 

SUBROUTINE  MATEQN3  (XI .  X2,  X3,  Y1 , Y2, Y3, A1 1 , A12, A13, A21 , A22, A23, 

*  A3 1 i A32 , A33 , LC , LN , LQ ) 

C 

c********************************** ****** ***************************** 

c 

C  THIS  SUBROUTINE  SOLVES  THE  THREE  SIMULTANEOUS  BAND  MATRIX 

C  EQUATIONS 

C 

C  All *X1  +  A12*X2  +  A13*X3  ■  Y1 

C  A21IX1  +  A22*X2  +  A23*X3  -  Y2 

C  A31IX1  +  A32*X2  +  A33*X3  -  Y3 

C 

C  FOR  XI,  X2,  AND  X3 

C 

C  A! J  ARE  9  BAND  MATRICES  OF  LENGTH  LQ,  WORKING  LENGTH  LN. 

C  AND  WIDTH  LC 

C  (THESE  MATRICES  ARE  ASSUMED  TO  BE  CORNER  ADJUSTED,  I.E.  THE 

C  CORNER  ELEMENTS  ARE  STORED  IN  (1.1)  AND  (LN, LC) .  ETC.) 

C 

C  XI  AND  Y I  ARE  VECTORS  OF  LENGTH  LQ  AND  WORKING  LENGTH  LN 
C 

c**************** ********************************************** ******* 

r 

DIMENSION 

*  XI (LQ) , X2(LQ) , X3(LQ) , Y1 (LQ) ,  Y2 (LQ) , Y3 (LQ) . 

*  All (LQ,LC) , A12(LQ,LC) , A13(LQ,LC) , 

*  A21 (LQ, LC)  A22 (LQ, LC) , A23 (LQ, LC) , 

*  A31 (LQ, LC) , A32 (LQ, LC) , A33 (LQ, LC) 

INITIALIZATION 


LP-LN+1 

L-(LC-l>/2 

LM*LN-L-1 

IF (LC.GE.LN)  L»LN 

DO  3  I“1 ,LN 

XKI)-Yl(I) 
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> 


I"  ■ 


X2  < I > -Y2 1 1 > 

X3<I)-Y3(I) 

3  CONTINUE 

C 

c  DOWNWARD  GAUSSIAN  ELIMINATION  WITH  PIVOTING 


DO  401  K'»i,LN 
IF (L.EQ. LM)  L«LN 
IF(L.LT.LN)  L"L+1 
C 

U-ABS<A11 <«,!)) 

I-K 

h«i 

DO  113  J*K,L 

IF(J.EQ.K)  GO  TO  111 

V-ABS(AU  (J,  1) ) 

IF(V.LE.ll)  GO  TO  111 

U«V 

M*1 

I*J 

111  V-ABS(A21<J,l)> 
IF(V.LE.U)  GO  TO  112 
U-V 

M*=2 

I«J 

112  V»ABS(A3l ( J,  1 ) ) 
IF(V.LE.U)  GO  TO  113 
u*v 

M=3 
I  ■*  J 

113  CONTINUE 
IF(I.EQ.K)  GO  TO  115 
IF(M.NE.l)  GO  TO  116 
DO  114  J«1,LC 

U*A1 1 <K, J) 

All (K, J)«A11 (I, J) 
A11<I,J)-U 
U*A12(K, J) 

A12<K, J)«A12(I, J) 

A12( I, J) *U 
U-A13<K,  J) 
A13(K,J)«A13(I,  J) 
A13(I,J)-U 

114  CONTINUE 
U«X1<K> 

Xl<IO-Xl(I) 

X1<I)-U 

GO  TO  120 

115  IF (M.EQ. 1)  50  TO  120 

116  IF <M. NE. 2)  GO  TO  1  IB 
DO  117  J«1,LC 


^ 


) 

'i 

\ 


i 

j 


t 
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U-AIKK.J) 

A11<K,J)-A21<I,J> 

A21(I,J)-U 

U-A12<K,J> 

A12(K,J)-A22(I, J) 

A22 ( I i J) »U 
U-A13(K,J) 

A13(Kf J)-A23(I, J) 

A23  < I , J  > -U 
117  continue; 

U-X1<K) 

XI <K> -X2  < I ) 

X2(I)-U 
GO  TO  120 

US  DO  119  J-1,LC 
U-All <K,J) 

Ail(K, J)-A31(I, J) 

A31<I.J)-U 
U-A12<K, J) 

A12(K,J)»A32(I, J) 

A32(I, J)-U 
U-A13<K,J> 

A13(K,J)-A33(I, J) 

A33 ( I , J) *U 

119  CONTINUE 
U-Xi <K) 

XI <  K  > *X3 ( I ) 

X3(I)-U 

120  CONTINUE 
C 

DO  126  I-K,L 
IF(I.EQ.K)  GO  TO  123 
U-A11<I,1>/AU(K,1) 

DO  122  J*1 , LC 

IF(J.NE.l)  A13 ( I, J-l J "All < I, J) -All <K, J) *U 
All (I, J)«A12<If J)-A12(Kf J)*U 
A12(I, J)-A13(I, J) -A13  <K, J)*U 

122  CONTINUE 
A13(I,LC>-0. 

XI < I ) -XI < I > -XI <K)*U 

123  CONTINUE 
U«A21(I,1)/AU(K,1) 

DO  125  J-1,LC 

IF(J.HE.i)  A23 ( I , J-i ) "A21 < I , J) -Ai 1 <K, J) *U 
A21(I,J)-A22(I, J)-A12<K,J>*U 
A22< I , J)«A23( I , J) -A13  <K, J) *U 
123  CONTINUE 

A23(I,LC)"0. 

X2(I)"X2(I)-X1(K)*U 
U"A31  <1, 1)  /All  <K,  1) 

DO  127  J*=l ,  LC 

IF ( J. NE. 1 )  A33 ( I , J— 1 > »A31 < I , J) -A1 1 (K, J) *U 
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A31 <ZP J>>A32(Zf J)-A12<KP J)*U 
A32(If J)-A33(I,J)-A13(Kf J)*U 

127  CONTINUE 
A33(I,LC>-0. 

X3 ( I  )*X3 ( I ) -XI (K) IU 

128  CONTINUE 
C 

U"ABS(A21 (K, 1 ) ) 

I-K 

M-2 

DO  213  J-K,L 
IF(J.EQ.K)  GO  TO  212 
V«AB8(A11 (J, 1 ) ) 

IF  <V. LEf  U)  GO  TO  211 

U«V 

M»1 

I«J 

211  V*>ABS (A21  < J, 1 ) ) 

IF (V( LE. U)  GO  TO  212 

U*V 

M=2 

I«J 

212  V*ABS(A31 (J, 1 ) ) 

IF(V.LE.U)  GO  TO  213 
U*V 

M=3 

I«J 

213  CONTINUE 
IF(I.EQ.K)  GO  TO  215 
IF(ri.NE.2)  Gu  TO  2it 
DO  214  J-l.LC 

UnA21 (K, J) 

A21 (K, J)«A21 (I, J) 

A21<I,J)-U 
U-A22(K, J) 

A22(K,J)=A22(I,J> 

A22 ( I , J )  »U 
U«A23<K, J) 

A23(K, J)*A23(I, J) 

A23U, J)«U 

214  CONTINUE 
U«X2 (K) 

X2(K)-X2(I) 

X2(I)-U 

GO  TO  220 

215  IF (M. NE. 3)  GO  TO  220 

216  IF (M. NE. 1 )  GO  TO  218 
DO  217  J-1,LC 
U«A21(K,J> 

A21 (K, J)-A11 (If J) 

All (I,J)-U 
U«A22(K, J) 


rr;v.  *  - v  ■ 


A22 (K,J)=A12(I,J) 

A12(I, J)-U 
U«A23(K, J) 

A23(K,J)-A13(I,J) 

A13( I, J)=U 

217  CONTINUE 
U«X2 (K) 

X2(K)-X1(I) 

Xl(I)-U 

GO  TO  220 

218  DO  219  J*1,LC 
U«A21(K,J) 

A21 <K, J) »A31 < I , J> 

A31 <1, J)-U 
U«A22(K,J) 

A22(K, J)«A32(I, J) 

A32(I,J)-U 
U-A23(K, J) 

A23(K, J)-A33(I, J) 

A33 <  I , J) *U 

219  CONTINUE 
U-X2(K> 

X2(K)«X3(I) 

X3(I'«U 

220  CONTINUE 

i 

DO  228  I-K.L 
IF(I.EQ.K)  GO  TO  223 
U*All<I.li/A21(K, 1) 

DO  222  J*15LC 

IFtJ.NE. 1)  A13(I,J-1)*A11 (I, J)-A21 (K, J)*U 
All  <1,  J)=A12(I,  J)-A22(K,  JXU 
A12(I,J)=A13(I,  J)-A23(K.  J)*U 

222  CONTINUE 
A13 ( I , LC) =0. 

XI < I > «X1 (I) — X2  <K> tU 
U»A21 (I, 1)/A21 (K, 1) 

DO  225  >1,LC 

IF (J.NE. 1)  A23 < I , J-l ) «A2i < I ( J > -A21 <K, J)*U 
A21 ( I , J ) «A22 <  I ,  J  > -A22 (K, J ) *U 
A22(I,  J )  -A23 < I  p  J)”A23(K,  J)  *U 
225  CONTINUE 

A23 ( I , LC) «0. 

X2(I)«X2(I)-X2(K)3iU 

223  CONTINUE 
U»A31(I,1)/A21(K,1> 

DO  227  J-1,LC 

IF  (J.NE.  1)  A33(I,  J-D-A31  (I,  J)-A2i  (K,  J)*U 
A31<If J)“A32(I, J)-A22(K,J)*U 
A32(I, J>-A33(I, J)-A23(K,J)$U 
227  CONTINUE 

A33 ( I , LC) =0. 


XS ( I ) -X3  ? I ) -X2 (K) tU 

228  CONTINUE 
C 

IF <K.EQ. LN)  GO  TO  401 
U-ABS<A31<Kf 1) ) 

I-K 

M-3 

JL-K+1 

DO  313  J«JL,L 

V»ABS(Ali(J,l)) 

IF(V.LE.U)  GO  TO  311 

U-V 

h«l 

I-J 

311  V-ABS(A21 (J, 1) ) 
IF(V.LE.U)  GO  TO  312 
U-V 

M-2 

I-J 

312  V»ABS(A31(J,1>) 
IF(V.LE.U)  GO  TO  313 
U-V 

M-3 

I-J 

313  CONTINUE 
IF(I.EQ.K)  GO  TO  320 
IF (M. NE. 3)  GO  TO  316 
DO  314  J-l.LC 

U—  A  *T  4  i  %J  H 

1  \^,U/ 

A31<K,J)-A3i(I,J> 

A31 (I, J)-U 
U-A32(K, J) 

A32(K, J)-A32<I,J) 

A32 ( I , J) -U 
U»A33<K* J) 

A33<K, J)-A33(I, J) 

A33  < I f  J) «U 
314  CONTINUE 
U-X3  <K) 

X3<K)-X3(I> 

X  3 ( I > -U 
GO  TO  320 

316  IF (M. NE. 1 >  GO  TO  318 
DO  317  J«1,LC 
U-A3KK,  J) 

A31  <K,J)-AU<I,  J) 

All (I, J)-U 
U»A32(K, J) 
A32<K,J)»A12<I,J) 
A12(I,J)-U 
U-A33<K, J) 

A33(K, J)«A13(I, J) 
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UPWARD  GAU88IAN  ELIMINATION 


L*1 

DO  507  K»1,LN 
I-LP-K 

U-X3(I> 

IF (I.EQ.LN)  GO  TO  502 
DO  501  J«2,  L 
IJ-I+J 

U-U-A32  <I,J-1)*X1<IJ-1) -A33 ( I , J-l > 1X2 ( I J-l > -A31 <1, J) *X3 ( I J-l ) 
IF(L.GE.LC)  U*U-A32 < I , LC) 1X1 < I+LC) -A33 ( I , LC) tX2  < I+LC) 

502  X3(I)-U/A31(I,1) 

U-X2 ( I ) -A22  < 1 , 1 ) >X3 ( I  > 

IF(I.EQ.LN)  GO  TO  504 
DO  503  J«2,L 
IJ=I+J 

503  U«U-A23<I,  J-l )  *X1 < I J-l ) -A21 (I ,  J)  *X2 ( I J-l > -A22  ( I , J) *X3(I J-l > 

I F ( L . GE . LC )  U-U-A23 ( I , LC ) *  X 1 < I +LC ) 

504  X2< I ) -U/A21 (I, 1) 

U-X1(I)-A12<I,1>*X2<I)-A13(I,  1>*X3<I) 

IF(I.EQ.LN)  GO  TO  506 
DO  505  J=2, L 
IJ=I+J 

505  U-U-All  (I,  J)*X1  <IJ-l)-A12<I.J)tX2(IJ-l>-A13<I,  J)*X3(IJ-1> 

506  X1(I>=L7A11<I,1) 

IF(L.LT.LC)  L-L+l 

507  CONTINUE 

RETURN 

END 

SUBROUTINE  RUFVAR (X, XNUE) 

COMMON  G,  PR,  REY,  XMINF,  OMEGA,  BO,  TW,  P10,  TIC,  RIO,  VIS10,  TE, 

1  PE,  RE, UE, VISINF, SU, EPS, DS, DYW,  SI , ERROR,  TC, TA, IEDGE, IEND1 ,  INTACT. 

2  PRT , XXK, BTRX, XLAM, VARPRT, XINTER,  8EP0, ICHS (8) , IPRN (9) , EO (300) . 

3  EN  (300) ,  EP  (300) ,  ETO  (300) ,  ETN  (300 ) ,  ETP  (300) ,  FO  (300) ,  FN  ( 300) ,  J2DA, 

4  FP ( 300 ) , TN ( 300 ) , TO  <  300 ) , XNN ( 300 ) , VN ( 300 ) , VO ( 300 ) , VP ( 300 ) ,  TP ( 300 ) , 

5  D1 (300) ,  D2 (300) , D3 (300) 

COMMON/BLRVAR/BYO (300) , BYN (300) ,  BYP (300) , XFY(300) , XSS(300) , 

1  XBRE, XHRE, XLS, XCD 
DIMENSION  XD(300) 

DATA  PI/3.1415926535/ 

mt*»mmmm*m*»«*m*****tmm«t 

INPUT  DATA  DIMENTIONALIZED  BY  L,  THE 
LENGTH  OF  THE  MODEL . 

mtMmmmmtmt  mm*  tit  *********** 

XSR-O. 


1  39 


v  v  siv 


C 


c 

c 


94 


95 


96 


XL8-.0046 

XLD8-.0046 

XBRE-.0023 

XDRE-. 0023 

XHRE-.00115 

XCD-0.6 

XD0-0. 00023 

NM1-IEDGE-1 
DO  99  N“I , NM1 
XBREE- XBRE 

XBREE  IS  A  DUMMY  VARIABLE 

XYY-XNN(N)*EPS 
IF(XSR.EQ.l)  GO  TO  95 
IF (XYY.LE. XHRE)  GO  TO  94 
XBREE-0.0 

BYP<N)«1.0-(XBREE*XDRE)/<XLS*XLDS> 

XFY <N> ■ ( 1 . O-XBREE/XLS) /BYP  <N) 

XD(N) -XBREE 
GO  TO  96 

DFUNCT=SQRT<  <XD0**2> /4. -XYY442) 

XD (N/ -DFUNCT 

BYP(N)«1.0-(PI*XD<N>**2) /(4.*XLS*»2) 

XFY  <N>  °  < 1 .  O-XD (N> /XLS) /BYP (N) 

IF  ( J2DA.  EG).  1 )  GO  TO  97 

XSS(N)=(XCD*XD<N>*X) / (BYP(N) *XLS**2*RE*UE*XNUE: 
GO  TO  99 


97  X0S<N)-<XCD*XD<N)*X)/<S*BYP<N)*XLS**2*RE*UE*XNJE> 
99  CONTINUE 
RETURN 
END 
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Appendix  F 


Four  Key  Subsystems  Within  Computer  Code 

Nondimensionalizing  the  Variables  and  Initializing  the  Grid 
Prior  to  entering  the  computational  loop  the  working 
variables  were  nondimensionalized  or  normalized.  These 
variables  were  listed  below  along  with  a  oo’-inition  of  each 
The  format  selected  was  to  present  the  coded  variable  on 
the  left  side  of  the  equal  sign  and  the  real  or  physical 
definition  on  the  right  side  of  the  equal  sign.  No  explan¬ 
ation  was  included  as  to  choice  of  normalizing  factors. 


T1  0 


1 


(y-dm; 


T 

o 


T  (y-1) 

oo  x  1 


,,2 

IS 


CO 


R1  0 


T  (y-1  )M‘ 

'  1  r. 


(F. 
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.  v.  ^  '.•;  .  s  . 


TINF  =  1 


T  /T 
o  00 


T  (y-I)M' 


(y~dm: 


With  Eq  (F.l)  defined  for  all  cases,  some  others  depended 


on  the  value  of  n 


Vli’j.0  = 


were  not  equal  to  zero,  then 


u 


[ (y-1  )m2  ] 


=  [  (y-i)m^  ]  w/2 


(Re  ) 


VISIKF  = 


[(Y-l)Kf  ]  '[Cf] 


where  the  reference  temperature  was  taken  as  T 

However,  for  the  case  where  w  was  equal  to  zero,  the  quan 
tities  of  Eq  (55)  plus  one  were  defined  as  follows: 


TC  =  S  =  1 98.6 
T  (y-l)M2  Tref 


VISIO  = 


.1  [ 


r:  (y-i)i-F 

m  x  1  '  n 


.  r  To  1  1-5r  Tref  +  198-61  .  % 
LTrafJ  L  ^  rT9°J  uref 


VISINF 


r 

■  (T  + 

V  OO 

1 98.6) 

(y-I)M^ 

1.5  -I 

1/2 

EPS  = 

(T  (y 

OO  ' 

-i  )m* 

OO 

+ 

198.6 

l 

- 

Re 

00 

- 

(F.3) 

rTref  1 

1.5 

’  T 

00 

+ 

198.61 

1/2 

1/2 

L  Tco  J 

T 

L  ref 

+ 

198.6J 

= 

f!Jref  /u«5 

Re 

00 

- 

L  J 

(F.  4 ) 

=r 

T 

OO 

v-5r 

T 

ref 

+ 

108.6  I 

[ 

T 

ref  . 

J  L 

f 

oo 

+ 

198.6  J 

These  quantities  wero  frequently  used  in  the  grid  computation 
and  provided  a  summary  of  the  nondimensional izing  techniques 
used  throughout  the  code.  Eefore  beginning  this  computation 
within  the  grid,  however,  there  had  to  be  an  initialization 
of  the  profile. 

Initialization  began  by  defining  Y  in  the  code  as  the 


distance  measured  along  the  n  axis.  Any  An,  .  was  defined 

J 

*  <| 

^  Aq •  which  yielded  a  fine  mesh  of  nodal  point 
0 


as  / AnK+1  V 

\  / 


near  the  surface  and  an  adequate  spacing  toward  the  edge. 

Y  values  were  assigned  by  successively  adding  all  An  values 
from  the  surface,  to  the  point  in  question.  Then,  three 
hypothetical  successive  columns  of  nodes  were  created  by 
the  following  statements: 
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D1  =  D2  =  D3  =  0.,  from  tho  surface  to  the  edge  of  the 
boundary  layer.  Incorporating  the  notation  of  fig  1, 


Y.  for  all  j  from  the  surface 

.] » 

to  the  edge  of  the  boundary 
layer . 

In  a  similar  manner,  three  successive  stations  of  F,  0_,  e, 
and  c  were  assigned  values  of  1.0.  Finally,  all  coeffi¬ 
cients  of  the  system  of  finite  difference  equations  were 
set  equal  to  0. 

This  initialization  provided  the  primer  to  begin  the 
backward  differencing  along  the  £  direction  and  the  central 
differencing  along  the  n  direction.  The  finite  differencing 
system  v/as  unconditionally  stable  for  increments  of  An  and 
A£,  and  the  iterative  stepping  procedure  along  £  damped  the 
error  due  to  the  grid  initialization  within  a  few  steps. 

Subroutine  Reystr 

This  routine  was  called  from  the  main  program  at  each 
station,  s^,  at  and  beyond  tho  point  of  transition  to  tur¬ 
bulence.  The  purpose  of  this  subroutine  was  to  calculate 
an  eddy  viscosity  for  the  inner  and  outer  regions  of  the 
two-layer  turbulent  boundary  layer  model. 

Computation  within  Reysir  began  with  Taylor  series 
expansions  of  F  to  the  third  order  partial  term  about  the 
first  station  at  the  wall.  With  values  for  F^_1f  F._9, 

and  F.  a  four-point  finite  difference  expression 


V.  .  =  V.  ,  =  V,  o 

i,j  1-1. j  1-2. j 


was  formed  for  _3F  and  the  coefficients  of  the  F  terms  at 

3p  w’ 

each  node,  one  through  four,  were  represented  hy  Y11,  Y12, 
Y13,  Y14  in  the  code.  Next,  a  nondimens ional  molecular 
viscosity-density  term  was  calculated  for  the  wall  with  a 
shear  stress  term  that  followed: 


XLM1W 


(F.  5 ) 


PI2 


9F 


w 


(F .  6 ) 


An  iterative  loop  was  begun  to  generate  the  nondimen- 
sionalized  inner  eddy  viscosity  model,  inner ,  of  Cebeci- 
Smith-Mosinskis  for  each  node  in  the  n  direction  for  the 
current  s^.  In  the  actual  code  and  the  following  the  calcu¬ 
lation  of  a  number  of  interim  quantities  that  did  not  neces¬ 
sarily  represent  any  real  boundary  layer  characteristic, 
three  important  computations  were  made.  First,  6/L  was 
calculated.  Next,  an  intermediate  quantity,  DD,  to  be  used 
later  in  the  outer  eddy  model,  was  calculated.  Finally, 

PI1 ,  another  intermediate  quantity  used  in  the  inner  model, 


was  computed: 

6/L  =  XNN . 

t) 


A  I"  —  ' 

_UeJ  j-(j-l) 
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edge  of  the  r 
boundary  layer 
DD  =  %  1 

j=2 


PI1  = 

2XRe 

OO 

edge  . 

£  Hui 

1  =  2  2 

[-T  .  T.n 

W  +  1  I 

T  T 

*-  e  e  J 

{(Y-OMV-5r  tjrwrr— - 

rii? 

LT»(Y-1)lF+198.6. 

uref  L  TeJ 

s 

where  the  s  was  included  for  the  case  of  conical  flow  only. 

Again,  a  _3F  term  was  generated,  but  using  only  a  three-point 

3n 

central  differencing  scheme  on  this  occasion.  The  final 

step  of  the  loop  was  the  actual  computation  of  e .  at 

^  r  _ inner 

the  current  node  j : 
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=  .  1  6 ( P1 1  )  (1  -  exp(-{  ( ? 1 1  )  (PI2)}  1/2/(26  ^2 


(pu) 


(F.6) 


t 


F(i, j+1  ) 


-  Yie  F(i,j) 


Y9 


F(i,  j-1  ) 


An , 


r (pp.k  i 


] 


where  Y8,  Y9»  and  Y10  were  coefficients  obtained  through 

Taylor  series  expansions  of  F(i,j-1)  and  F(i,j+1)  about  a 

point  F(i,i).  As  the  calculatio?i  of  c.  progressed  from 

the  wall  out  into  the  field  of  flow,  e.  ^  retained 

inner... 

_ 111 

U 

its  own  commuted  value  or  that  cf  r.  whichever  was 

inner. 

greater.  p 

The  outer  law,  £ou^er  was  computed  through  an  iterative 


U 


loop  similar  to  that  of  the  inner  model.  It  culminated 
with  the  expression 
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where  the  s  was  included  only  for  the  case  of  conical  flow. 

In  order  that  a  compatible  combination  of  computer  viscosities 
were  retained,  the  values  of  eddy  viscosity  from  the  outer  law 
replaced  those  of  the  inner  law  from  the  point  of  intersection 
of  the  graphs  to  the  edge  of  the  boundary  layer.  Graphically, 
this  was  depicted  in  Fig  3. 

Having  calculated  the  initial  eddy  values  for  the  inner 
and  outer  viscous  regions  of  the  boundary  layer,  it  was  appro¬ 
priate  to  subject  this  model  to  two  more  factors.  Both  were 
factors  of  degradation  and  were  included  to  better  describe 
the  character  of  turbulent  activity  within  the  boundary 
layer. 

Objections  have  been  raised  to  the  use  of  an  eddy  vis¬ 
cosity  term,  e,  in  place  of,  or  in  addition  to  the  molecular 
viscosity,  u,  of  a  fluid,  y  is  a  real  property  of  a  fluid, 
e  is  only  an  effective  description  when  a  fluid  is  in  motion, 
and  it  is  clearly  not  a  property  of  the  fluid.  But,  with 
reservation,  it  has  been  used  to  express  the  behavior  of 
turbulent  stresses  in  terms  of  mean  velocity  gradients  of 
a  flowing  fluid.  It  has  been  possible  to  obtain  a  satis¬ 
factory  description  of  mean  properties  within  turbulent 


flows  by  assuming  this  flow  to  behave  as  a  Newtonian  fluid, 
incorporating  an  eddy  viscosity  model  along  with  u»  and  in¬ 
cluding  two  factors  of  intermitten cy  when  appropriate  (Ref. 
38:25-26).  A  laminar  and  irrotational  flow  became  turbu¬ 
lent  as  it  passed  through  a  region  of  transition  in  which 
only  a  fraction  of  the  time  was  spent  in  a  turbulent  state. 
During  that  time  in  laminar  motion,  the  Reynolds  stress, 
hence  e,  would  have  been  zero.  Then,  to  adequately  describe 
the  effects  of  e  at  any  point  by  the  relative  fraction  of 
time  that  that  point  would  be  engulfed  in  turbulent  flow 
(Ref  34*117).  Therefore,  the  first  multiplicative  factor, 
called  an  intermittency  factor,  was  applied  to  e  to  more 
accurately  describe  the  £  within  the  transition  region. 

The  intermittency  or  probability  factor  cf  Dhsvar  and 
Narasimha  was  used  for  this  program.  The  factor  was 
computed  as  follows  (Ref.  13*28-29): 
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The  second  factor  was  then  considered.  It  was  observed 
by  Klebanoff  that  in  a  turbulent  boundary  layer  v;ith  a  free 
boundary,  as  the  free  stream  was  approached  the  turbulence 
became  intermittent.  This  intermittent  nature  was  observed 
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first  at  y/6  greater  than  .4  with  less  turbulent  intensity 
as  y/6  grew  larger.  It  was  thought  that  a  good  prediction 
of  turbulent  intensity  probably  depended  on  a  correct 
weighting  of  the  probability  density  for  the  turbulence 
of  the  free  stream  with  that  within  the  boundary.  It  was 
found  that  a  good  description  of  y*  was  a  Gaussian  integral 
curve  given  by 

Y’  =  1  (1-erfU'))  (F.-I2) 

where 

c'  '(Fi  y  (I  - 

These  expressions  indicated  that  the  edge  of  the  boundary 
layer  had  a  random  character  with  a  mean  position  at  y/o 
equal  to  .78.  The  edge  vacillated  from  y/6  equal  to  ,4 
to  y/6  equal  to  1.2.  Finally,  if  it  were  assumed  that  the 
free  stream  contributed  little  to  the  measured  turbulent 
quantities  of  the  boundary  layer,  an  allowance  could  be 
made  for  the  effect  of  intermitt ency  by  dividing  by  the 
factor  y*  (Ref.  28:15-19). 

Cebeci  used  the  approximate  expression  for  Eq  (F.12) 
to  give  a  multiplicative  version: 

/  -1 

y '  =  1  +  5.5  |  (Ref.  20:96)  (F.14) 

which  led  to  the  coding  for  this  second  factor.  If  Y1  were 
not  included,  then  a  newly  defined  viscosity  was 


78 J  -  5^ 


(F.13) 
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(F.15) 


‘  -  i  ♦  *  r(s) 

Including  Y'»  Shang  formed  the  following  model: 

f  u  1  (F.16) 

L2.75J 

For  purposes  of  this  study  Eq  (F.15)  became  eddy  model  zero, 
and  Eq  (F.16)  became  eddy  model  one.  Then,  whether  or  not 
Y’  was  included,  the  quantity  e  was  defined  by 

6  =  1  +  — -  (e-1)  (F.17) 


e  =  1  + 


1.75 


+ 1 


In  a  final  note,  the  decision  of  whether  to  use  eddy  model 
zero  or  eddy  model  one  depended  on  the  original  assumption 
that  either  the  free  stream  turbulence  had  an  effect  on  the 
e  of  the  boundary  layer,  or  it  did  not.  This  factor,  y', 
was  to  have  a  definite  effect  on  the  analytical  results, 
and  this  entire  subroutine  was  included  with  program  listing 
of  Appendix  E. 


*\ 1 


Subroutine  Cfstr  i 

Like  Reystr  this  routine  was  called  from  the  main  pro¬ 
gram,  But  unlike  Reystr,  Cfstno  performed  its  computation 
throughout  the  laminar,  transition,  and  turbulent  regions 
of  flow.  The  purpose  of  this  routine  was  tc  calculate  a 
Stanton  number,  a  measure  of  heat  transfer;  the  local 
coefficient  of  friction,  indicative  of  shear  stress  at  the 
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surface;  and  Reynolds  numbers  based  on  displacement  thick- 
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ness  and  momentum  thickness. 


(pu), 


Computation  began  with  coded  XLM1  in  the  program 

T^re 


The  formula  by  which  XLM1  was  computed  depended  on  the  value 


of  the  exponent  in  the  viscosity  law  of  Sutherland,  the 


value  of  this  exponent  being  specified  by  the  programmer, 


If  the  exponent  were  zero,  then 


XLM1  = 


r  h  ?/2rvi98-6i 

L  Te  J  LV198-6J 


If  this  exponent  were  one,  then  XLM1  was  one.  Otherwise, 


vtm-1  T 

XLM1  =  w 


( F .  1 


Next  to  be  calculated  were  'ransf ormed  quantities  similar  to 


q  or  heat  flux  and  r  or  shear  stress.  First,  the  same  four- 


point  finite  difference  scheme  used  in  Reystr  for  JF  was 

an  w 


repeated  at  this  point  to  calculate  3F  and  56_  Then 

9q  w  3n  w* 


the  transformed  t,  coded  TAU,  was  computed: 
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where 


TAUR  = 


«»-* "*  ,F- 


r\v  is  the  value  of  n  at  v  =  h  (height  of  the  roughness 


element )  . 


The  skin-friction  coefficient  is  increased  by  the 
addition  of  the  roughness  element  and  it  may  be  defined  as: 


T  +D/.BC 
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fw - 5— 
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where  D  =  /b  1  u2  C-.  dy  dz.  x  and  D/BC  may  be 

o  o  ~  D  J  J 


sw 


written  in  dimensional  variables  as  follows 


t  +  D  =  u  3  u 
w  BC  w  3  y !  w  2  BC  0 


+  1  CDb  fk  2  , 

/  u  dy 


The  dimension  of  the  element  is  constant  along  y  for 
the  case  of  rectangular  roughness  elements.  Non-dimension- 
alizing  the  variables  and  stretching  the  y  coordinate. 
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Converting  to  transformed  coordinate  n , 
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In  the  program,  the  increased  t  cue  to  roughness  elements 
is  defined  as  TAUR ,  by  equation  (?.21). 

Following  t,  the  transformed  q,  coded  QS,  was  replaced  by 
the  following  expression: 


QS  = 
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For  the  case  of  the  axisymmetric  flow,  both  TAU  and  QS  were 
derived  by  the  non-dimensional  station,  s^.  With  this, 
preliminary  calculations  were  completed. 

A  Stanton  number  and  coefficient  of  friction  followed 
next  in  the  computation.  If  equaled  T  ,  there  was  no 
heat  transfer  and  ST,  coded  STMO,  was  zero.  Otherwise, 


The  model  from  which  this  expression  came  was 


Ste  =  q 

pP«err:a-h,.j 


(F.28) 


For  the  calculation  of  of|local  station.  ^oded  CFNO, 


CFNO  =  p  Uref  ‘  P°°U°°L  ^pv\i  /Ue\2  (2X)'^  3J: 
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With  St  and  c„  computed,  only  the  transformed  expressions 

Ilocal 

for  Regtt  and  ReQ  remained.  Coded  as  REYDT  and  REYMT ,  these 
quantities  were  computed  from  the  following  statements: 


REYDT 


=  peuexreal  6* 
L 


REYMT  =  peuexreal  0 


(F.3 


Subroutine  RUFVAR 

The  effect  of  surface  roughness  on  compressible  turbu¬ 
lent  boundary-layer  is  modelled  by  distributed  sources  and 
sinks  and  blockage-terms  in  the  appropriate  governing 
equations.  This  subroutine  was  called  from  the  main  pro¬ 
gram  at  each  station,  s.,  throughout  the  laminar,  transition 
and  fully  turbulent  regions  of  the  flow.  The  purpose  of 

this  subroutine  was  to  calculate  the  source/sink  term  ($  ) 

s  s 

and  the  blockage  terms  (f(y)  and  B(y)). 

To  initiate  this  subroutine,  the  following  quantities 
were  specified  as  input. 

XSR  -  a  flagged  ouantity  to  specify  whether  the 
roughness  elements  are  of  spherical  or 
rectangular  shape 

XBRE  -  the  breadth  of  the  rectangular  element 

XDRE  -  the  depth  of  the  rectangular  element  in  the 
direction  of  flow 

XHRE  -  height  of  the  roughness  element 

XLS  -  centre  to  centre  spacing  between  adjacent 
elements  facing  the  flow 

XLDS  -  centre  to  centre  spacing  between  adjacent 
elements  in  the  direction  of  the  flow 


XCD  -  coefficient  of  drag 

XDO  -  the  diameter  of  the  roughness  element  with  circular 
cross  section  at  y  =  0 

DFUNGT  -  the  shape  of  the  sph  erical  element 

The  detailed  dimensions  of  the  roughness  elements  for 
rectangular  cross-section  are  given  in  Fig  5.  The  input 
data  is  non-dimensionalized  by  the  length  of  the  model. 

After  the  size,  shape,  and  the  spacing  of  the  roughness 
elements  are  specified,  the  subroutine  calculates  the  block¬ 
age  terms  f(y)  and  B(y)  along  p  as  follows: 

(a )  For  Elements  with  Rectangular  Cross-Sections 


BYP(N)  =  1.0  -  XBRE  *  XDRE 

XLS  *  XLDS 

and 


XFY(N)  =  1.0  -  ( XBRE ) / ( XLS ) 

BYP(N) 


( b)  For  Elements  With  Circular  Cross-Sections 

BYP(N)  =  1.0-Ti  XD(N)2 

4XLS2 

and 

XFY(N)  =  1.0  -  XD(N)/XLS 

BYP(N) 


After  calculating  ohe  blockage  terms,  the  subroutine 
further  computes  the  source/sink  term  as  follows: 


XSS(N)  = 


For  elements  with  rectangular  cross-section,  D(y)  was 
set  equal  to  XBFE. 

This  completed  calculations  within  this  routine,  and 
further,  completed  the  formal  description  of  four  impor¬ 
tant  subsystems  within  ITRACT.  Again,  this  subroutine 
was  included  with  the  program  listing  of  Appendix  E.  In 
this  appendix  consideration  was  given  to  the  important 
concepts  of  the  nondimensionalizat ion  of  working  quanti¬ 
ties,  and  initialization  of  the  grid.  Also  included  was 
a  brief  description  of  the  three  subroutines  used  in  the 
computation  of  eddy  viscosity,  heat  transfer,  and  skin 
friction  and  roughness  variables.  The  theory  presented 
in  this  appendix  should  provide  a  center  understanding  o 
the  code  in  general,  and  the  modification  for  surface 
roughness  specifically. 
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to  account  for  the  blockage  effect  due  bo  the  roughness  elements.  The 
modified  governing  equations  were  then  transformed  using  the  transformation 
of  Probstein-Elliott  and  Levy-Lees.  The  resulting  equations,  with  appro¬ 
priate  boundary  conditions,  were  solved  by  finite-difference  techniques  to 
determine  the  non-dimensional  velocity  components  and  temperature  at  a 
finite  number  of  nodes  in  the  boundary-layer  field  of  the  flow. 

To  establish  the  authenticity  of  the  original  code,  ITRACT,  some  smooth 
surface  results  were  first  computed  and  compared  with  the  experimental  data 
of  Dr.  Fiore  and  Dr.  Cole  for  turbulent  flows  over  smooth  surfaces.  After 
the  accuracy  of  the  code  was  established  for  smooth-surface  calculations, 
the  code  was  modified  in  order  to  render  it  capable  of  predicting  the 
influence  of  surface  roughness  on  compressible  turbulent  flow.  The  modified 
code  was  then  used  to  obtain  results  for  rough-surface  boundary  layers  and 
Pj|he  computed  results  were  compared  with  the  experimental  data  of  Dr.  Fiore 
'■i’<or  the  case  of  supersonic  flow  over  a  rough  flat  plate.  The  agreement 
between  the  computed  and  the  measured  velocity  profiles  was  quite  satis¬ 
factory.  The  corresponding  temperature  profiles  agreed  well  everywhere 
except  very  near  the  wall;  a  possible  reason  for  this ■ discrepancy  is  offered 
in  this  study.  Unlike  the  previous  studies  of  rough-surface  boundary  layers, 
the  present  study  makes  no  modification  to  the  ‘turbulence  model  employed. 
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